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CHAPTER  1 


INTRODUCTION 


1.1  Limit  Cycle  Oscillations 

Limit  Cycle  Oscillations  (LCO)  are  known  to  occur  on  various  operational  vehicles.  In 
particular,  LCO  has  been  a  prevalent  aeroelastic  problem  on  several  current  fighter  aircraft.  It 
usually  occurs  for  aircraft  with  external  stores  throughout  the  transonic  flight  regime  (Refs  1-4). 
Complicated  by  the  aircraft-store  system,  its  mechanism  still  remains  to  be  fully  understood. 
Meanwhile,  a  business  jet  wing  LCO  was  also  reported  recently  (Refs  5,6).  This  has  been  a 
serious  concern  since  there  are  few  analytical  techniques  available  for  LCO  prediction  and  an 
insufficient  understanding  of  its  physics. 

LCO  can  be  characterized  as  sustained  periodic  oscillations  which  neither  increase  or  decrease  in 
amplitude  over  time  for  a  given  flight  condition.  Using  an  s-domain  unsteady  aerodynamic 
model  of  the  aircraft  and  stores,  Chen,  Sarhaddi  and  Liu  (Ref  7)  have  shown  that  wing/store 
LCO  can  be  a  post-flutter  phenomenon  whenever  the  flutter  mode  contains  low  unstable 
damping.  This  type  of  flutter  mode  is  called  a  “ hump  mode”.  Since  the  aircraft  structure  usually 
contains  structural  nonlinearity,  such  as  friction  damping,  this  amplitude-dependent  friction 
damping  can  suppress  the  growth  of  amplitude,  thus  resulting  in  a  steady  state  oscillation.  This 
is  known  as  the  nonlinear  structural  damping  (NSD)  model  of  the  wing/store  LCO.  Although  not 
thoroughly  proven  through  tests  or  numerical  simulations,  results  of  the  NSD  show  excellent 
correlation  with  flight  test  LCO  data  of  F-16  throughout  subsonic  and  transonic  Mach  numbers. 

On  the  other  hand,  other  researchers,  notably  Cunningham  and  Meijer  (Ref  8)  believe  that  the 
wing/store  LCO  is  due  largely  to  transonic  shock  oscillation  and  shock  induced  flow  separation. 
This  is  called  the  Transonic  Shock/Separation  (TSS)  model.  Edwards  has  suggested  the  TSS 
model  and  viscous  effects  are  two  major  factors  that  cause  transonic  LCO  for  wings.  He  also  has 
studied  the  shock  buffet  phenomenon  in  addition  to  transonic  LCO  (Ref  9).  It  should  be  noted, 
however,  that  there  is  no  conflict  between  the  NSD  model  and  the  TSS  model  in  that  both 
physical  effects  may  contribute  to  LCO. 

Note  that  the  method  in  Ref  7  used  a  damping  correlation  technique  for  LCO-onset  prediction, 
whereas  that  of  Ref  8  is  a  semi-empirical  method  that  requires  steady  and  unsteady  data  input 
from  wind  tunnel  measurement. 

1.2  LCO  Prediction  Methods 

Recent  renewed  interest  in  LCO  is  motivated  by  the  need  to  better  predict  and  understand  fighter 
LCO  and  also  by  the  rapid  advent  of  CFD  methodology  in  aeroelasticity.  Currently,  there  are 
two  potential  computational  approaches  for  LCO  prediction/investigation:  the  time-marching 
approach  using  a  high-level  CFD  code  such  as  CFL3D  (Ref  10)  developed  and  supported  by 
NASA/Langley,  and  the  Frequency-Domain  POD/ROM  EigenMode  approach  (Ref  11) 
originated  by  Dowell  and  Hall  of  Duke  University.  The  former  is  a  conventional  time-domain 
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CFD  method  whereas  the  latter  is  a  frequency-domain  CFD  method  using  aerodynamic 
eigenmodes. 

For  nearly  twenty  years,  the  aerospace  industry  has  been  lacking  of  a  viable,  efficient  CFD 
method  for  transonic  flutter  applications.  The  time-marching  CFD  unfortunately  remain 
inefficient  as  a  practical  tool. 

Under  an  AFOSR/STTR  Phase  I  contract,  ZONA  has  recently  conducted  a  feasibility  study  on 
the  simulation  of  transonic  LCO  for  a  supercritical  NLR7301  airfoil  using  the  CFL3D/Navier- 
Stokes  code.  It  was  found  that  the  time-marching  LCO  solution  is  turbulence-model  dependent 
and  time-step  sensitive.  One  LCO  solution  case  typically  would  take  96  hours  on  a  1GHz  CPU 
computer. 

Earlier  work  of  Bendiksen  (Ref  12)  and  recent  work  of  Sheta  et  al  (Ref  13)  also  used  time¬ 
marching  CFD  approach  for  Euler  transonic  and  N-S  incompressible  LCO  simulations, 
respectively.  Based  on  the  cases  studied  by  ZONA,  the  time-marching  Euler  simulation  of  a 
transonic  LCO  requires  one  half  of  the  computing  time  required  of  a  Navier-Stokes  simulation, 
due  to  the  reduction  in  grid  points.  The  computing  time  required  for  incompressible  Navier- 
Stokes  LCO  simulation  is  about  two  orders  of  magnitude  less  than  a  transonic  Navier-Stokes 
LCO  simulation,  due  to  the  absence  of  transonic  nonlinearity,  thus  allowing  much  coarser  time 
steps  and  a  reduction  in  subiterations.  For  transonic  Navier-Stokes  simulation  of  LCO,  the 
requirements  in  time  steps,  grid,  subiterations,  etc.,  become  much  more  stringent,  resulting  in 
days  of  computing  time  for  one  LCO  case.  This  leads  to  the  conclusion  that,  if  a  transonic 
Navier-Stokes  simulation  were  demanded,  a  time-marching  approach  would  be  computationally 
inefficient  as  a  viable  approach  for  realistic  3D  transonic  LCO  prediction/investigations. 

On  the  other  hand,  the  Frequency-Domain  POD/ROM  (Proper  Orthogonal  Decomposition  / 
Reduced  Order  Model)  method  of  Duke  University  is  considered  a  recent  breakthrough  in  the 
transonic  computational  aeroelasticity.  Within  the  last  five  years,  the  Frequency-Domain 
POD/ROM  method  of  Dowell  and  Hall  has  been  proven  to  be  highly  efficient  and  shown  to  be 
widely  applicable  to  aeroelastic  problems  for  turbomachinery  cascades,  control  surface  freeplay, 
airfoil/wing  flutter  and  LCO  (e.g.,  Refs  11,  14-16,  Hall  and  Dowell,  3  Papers  SDM2001 
presented  in  Appendix  A-C  of  this  report).  As  evidenced  by  their  recent  studies,  under 
AFOSR/STTR  contract  support,  the  Frequency-Domain  method  with  the  Harmonic  Balance 
scheme  is  about  two-orders  of  magnitude  faster  than  conventional  time-marching  methods. 

For  example,  the  3D  Time  Linearized  code  can  generate  the  needed  aerodynamic  model  in  one 
day  for  a  given  Mach  number  and  flutter  points  take  less  than  1  minute  per  Mach  number  on  a 
workstation  computer. 

1.3  Frequency-Domain  POD/ROM  EigenMode  Approach 

Reduced-Order  Models  (ROM)  EigenMode 

The  eigenmodes  of  a  time-linearized  system,  which  may  be  thought  of  as  aerodynamic  states, 
were  computed  and  subsequently  used  to  construct  computationally  efficient,  reduced  order 
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models  of  the  unsteady  flow  field.  An  important  advantage  of  the  eigenmode  based  reduced 
order  modeling  technique  is  that  once  the  eigenmode  information  has  been  computed,  reduced 
order  models  can  be  constructed  and  used  to  calculate  the  response  at  different  frequencies  and 
mode  shapes  for  almost  no  additional  computational  effort.  Furthermore,  only  a  veiy  few 
eigenmodes  or  states  (degrees  of  freedom)  need  to  be  retained  in  the  model  to  accurately  predict 
the  unsteady  aerodynamic  response,  making  the  method  ideally  suited  for  rapid  flutter 
calculation  and  active  control. 

The  POD  Technigue 

The  Proper  Orthogonal  Decomposition  (POD)  technique,  also  known  as  Karhunen-Loeve  (Ref 
12)  expansions,  was  originally  introduced  to  determine  and  model  coherent  structures  in 
turbulent  flow  fields.  Using  the  POD  approach,  one  examines  a  series  of  “snapshots”  of 
experimental  or  computational  data,  each  at  a  different  instant  in  time.  These  solution  snapshots 
are  used  to  form  a  small  and  more  compact  eigenvalue  problem  that  is  solved  to  determine  a  set 
of  optimal  basis  functions  for  representing  the  flow  field. 

Frequency-Domain  POD/ROM 

Hall,  Thomas,  and  Dowell  (Ref  1 1)  developed  a  frequency-domain  form  of  the  POD  technique, 
and  applied  it  to  transonic  flows  about  airfoils.  In  particular,  they  used  a  time-linearized  CFD 
analysis  to  compute  unsteady  small-disturbance  flow  solutions  for  vibrating  airfoils  in  the 
frequency  domain  over  a  range  of  frequencies.  Basis  vectors  were  then  extracted  from  this 
frequency-domain  data  set  using  the  POD  technique.  The  resulting  basis  vectors  were  then  used 
to  construct  low  degree  of  freedom  reduced-order  models  of  the  unsteady  flow.  Finally,  the 
reduced-order  aerodynamic  model  was  combined  with  a  structural  dynamic  model  resulting  in  a 
compact,  but  accurate,  flutter  model. 

Harmonic-Balance  Method  for  Nonlinear  LCO  (Refs  14,  15) 

In  conjunction  with  the  frequency-domain  POD/ROM  method,  the  harmonic  balance  (HB) 
technique  can  be  used  in  the  treatment  of  nonlinear  LCO  problems  for  oscillating  wings  (with 
stores).  The  harmonic  balance  technique  is  used  to  recast  the  proposed  Euler  equations  into  a  set 
of  “higher-order”  equations  of  like  harmonics.  After  introducing  a  pseudo-time  term  into  the 
harmonic-balance  equations  so  that  they  may  be  solved  by  time  marching,  these  equations  are 
solved  by  conventional  CFD  techniques.  The  method  has  a  number  of  advantages  over  the  more 
conventional  time  domain  solutions.  Because  the  solutions  are  computed  in  the  frequency 
domain,  the  time-marching  algorithm  is  only  used  to  converge  the  solution  to  steady  state.  Thus, 
accelerating  techniques,  including  pseudo  time-marching  with  multi-grid  acceleration  can  be 
used.  Moreover,  for  problems  where  “engineering  accuracy”  is  required,  the  harmonic  balance 
series  can  be  truncated  into  just  a  few  harmonics.  The  result  is  that  the  proposed  harmonic 
balance  method  is  potentially  two-orders  of  magnitude  faster  than  conventional  time-marching 
methods  for  determining  LCO. 
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1.4  Merits  of  the  Frequency-Domain  Approach 

There  are  several  definite  advantages  of  using  the  frequency-domain  approach.  First,  the 
frequency-domain  harmonic  balance  method,  when  applied  to  high-level  equations,  can  be  at 
least  two-orders  of  magnitude  faster  than  the  nonlinear  time-domain  CFD  simulations.  Also,  this 
method  retains  essential  nonlinear  features  in  an  aeroelastic  system,  including  nonlinear 
structural  stiffness  and  damping  as  well  as  large  transonic  shock  excursions  including  viscosity. 
Second,  current  transonic  flutter  methods  using  time-domain/time-marching  CFD  in  conjunction 
with  various  versions  of  the  indicial  approach  (Batina,  Silva,  Refs  17,  18)  are  very  tedious  and 
computationally  expensive  and  their  accuracy  usually  depends  on  the  indicial  motion  imposed. 
By  contrast,  the  frequency-domain  formulation  of  POD/ROM  directly  solves  for  convenient 
aerodynamic  mode  shapes  which  can  be  stored  and  repeatedly  applied  for  a  large  range  of 
frequencies;  hence,  its  a  much  more  efficient  and  accurate  method.  Third,  the  frequency-domain 
approach  with  an  eigenvalue  solution  method  is  a  familiar  practice  to  the  structures/loads  and 
flutter  engineers.  The  proposed  method  after  ZONA’s  further  modification  is  expected  to  be 
well  accepted  within  the  industrial  environment. 

1.5  ZONA  Technology’s  R&D  in  LCO 

In  recent  years,  ZONA  has  been  extensively  engaged  in  the  R&D  of  LCO  including  its 
prediction  methodology,  control  and  F-16  and  F-18  LCO  data  correlation.  ZONA  has 
collaborated  closely  with  Lockheed-Martin/LMTAS  and  the  Seek  Eagle  office  of  Eglin  AFB  on 
F-16  wing/store  LCO  investigations.  The  results  including  the  definition  of  a  NSD  model  were 
presented  in  two  AIAA/SDM  papers  (Refs  7,19,  Chen  et  al,  Mignolet  et  al).  Under  a  recent 
NAVAIR  contract,  ZONA  also  worked  closely  with  team  member  Boeing/  St  Louis  on  a 
successful  R&D  in  the  reconfigurable  adaptive  control  of  the  F-18  LCO.  (Ref  20,  Nam  et  al). 
Supported  by  NASA/Langley,  onging  development  of  the  CFD/CSD  interfacing  using  BEM 
solver  (Refs  21,22,  Chen  et  al)  will  provide  a  new  spline  methodology  for  tightly  coupled 
aerodynamic-structural  interaction  for  3D  multiple  mode  LCO  studies  in  Phase  II. 

1.6  STTR  Phase  I 

In  August  1999,  the  AFOSR  awarded  an  STTR  Phase  I  contract  to  the  ZONA  Technology 
(ZONA)  team  to  develop  innovative  computational  aeroelasticity  methodologies  for 
understanding,  predicting  and  controlling  critical  nonlinear  aerostructural  interaction  (e.g., 
LCO/fl utter)  phenomena.  The  ZONA  team,  consisting  of  ZONA  Technology,  Inc.  and  Duke 
University  (Professors  Earl  Dowell  and  Kenneth  Hall),  has  successfully  accomplished  the  STTR 
Phase  I  contract.  The  recent  work  of  Duke  University  deals  with  the  further  generalization  of 
POD/ROM  EigenMode  method  for  2-D/3-D  flutter  and  LCO  solutions  (Refs  14,15,16). 

Under  this  STTR  contract,  ZONA  has  generalized  the  structural  compatible  model,  provided  a 
3D  flutter  solution  methodology  and  conducted  a  2D  CFD  simulation  of  transonic  LCO  (Ref  23). 


4 


1.7  Phase  I  Achievements 


Significant  milestones  have  been  reached  through  the  work  of  Phase  I.  These  are  in  four  related 
categories.  Firstly,  the  modeling  and  prediction  of  flutter  and  LCO  of  an  airfoil  and  control 
surface  with  structural  ffeeplay  in  2D  transonic  flow  has  been  accomplished.  Secondly,  the 
modeling  and  prediction  of  flutter  and  LCO  of  an  airfoil  undergoing  single  degree  of  freedom 
pitch  motions  due  to  nonlinear  aerodynamic  effects  associated  with  large  (inviscid)  shock 
motions  in  2D  flow  has  been  achieved.  Thirdly,  the  modeling  and  prediction  of  flutter  of  a  3D 
flow  about  an  elastic  wing  has  been  accomplished.  Note  that  in  the  latter  category,  we  are 
modeling  the  nonlinear  static  or  steady  flow  effects  in  the  transonic  regime  even  though  the 
aerodynamic  model  is  dynamically  linear,  i.e.,  the  shock  motions  are  assumed  to  be  proportional 
to  the  airfoil  motion  in  this  model.  And  finally,  a  2D  numerical  simulation  of  the  NLR  7301 
supercritical  airfoil  limit  cycle  oscillations  was  performed  using  the  time-marching  procedure  of 
the  CFL3D  code.  This  fourth  category  of  work  is  to  establish  benchmarks  for  the  viscous 
POD/ROM  methodology  developed  in  Phase  I  and  to  be  developed  in  Phase  II. 

In  Chapters  2-5  we  provide  a  summary  of  the  approach  and  results  obtained  in  each  of  these 
categories.  These  accomplishments  provide  a  firm  foundation  on  which  to  build  for  the  Phase  II 
work. 
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CHAPTER  2 


TRANSONIC  LIMIT  CYCLE  OSCILLATION  ANALYSIS  OF  AN 
AIRFOIL  WITH  CONTROL  SURFACE  FREEPLAY 


Summary 

The  transonic  flutter  and  limit  cycle  oscillations  of  an  airfoil  with  control  surface  freeplay  have 
been  determined  using  a  new  aerodynamic  modeling  technique  that  provides  greater  physical 
insight  and  understanding  by  tracing  the  true  root  locus  of  the  corresponding  linear  aeroelastic 
system.  This  in  turn  enables  a  very  computationally  efficient  harmonic  balance  technique  to  be 
used  in  determining  the  nonlinear  limit  cycle  oscillations. 

New  physical  insights  gained  include  the  rapid  change  in  flutter  mode  that  occurs  in  the 
transonic  Mach  number  range.  This  phenomenon  has  been  observed  in  experiments,  but  has  not 
been  previously  predicted  theoretically.  With  respect  to  LCO,  these  are  the  first  results  available 
in  the  transonic  range  for  the  configuration  studied.  The  model  also  predicts  significant  changes 
in  LCO  behavior  as  a  function  of  Mach  number.  However,  these  are  as  yet  unconfirmed  by 
experiments.  Up  to  high  subsonic  Mach  numbers,  the  flutter  and  LCO  results  are  similar  to 
those  previously  found  at  low  Mach  numbers. 


Introduction 


In  this  work,  we  will  assume  the  shock  motion  is  sufficiently  small  such  that  it  is  (linearly) 
proportional  to  the  airfoil  motion,  e.g.,  airfoil  motions  are  less  than  the  equivalent  of  one  degree 
in  angle  of  attack.  Using  an  Euler/CFD-based  reduced  order  aerodynamic  model,  a  thorough 
study  of  the  flutter  boundary  with  Mach  number  (M)  is  first  presented  in  the  absence  of  freeplay. 
Particularly  noteworthy  are  the  rapid  changes  of  flutter  modal  content  in  the  transonic  range. 
This  is  attributed  in  part  to  the  rapid  changes  of  center  of  pressure  location  as  the  mean  shock 
position  changes  with  Mach  number.  These  changes  in  the  modal  response  content  are  also 
found  in  the  limit  cycle  oscillations  (LCO)  which  are  encountered  when  control  surface  freeplay 
is  present.  Indeed  for  LCO,  the  modal  content  may  change  at  a  fixed  Mach  number  when  the 
dynamic  pressure  or  flow  density  is  varied. 

Below  M=0.80,  the  LCO  and  flutter  oscillations  are  qualitatively  similar  to  those  previously 
found  at  low  Mach  number  where  earlier  analyses  and  experiments  have  been  carried  out. 
However,  the  response  behavior  in  the  transonic  flow  regime  is  notably  different.  Of  special 
interest  is  the  occurrence  of  flutter  in  a  narrow  range  of  Mach  number  for  pitch  and  flap  (control 
surface)  dominated  motions.  Moreover,  beyond  a  certain  high  transonic  Mach  number  (after  the 
mean  shock  position  reaches  the  trailing  edge  of  the  airfoil),  neither  flutter  nor  limit  cycle 
oscillation  occurs. 
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Significance  ofLCO 


LCO  is  known  to  occur  on  various  operational  aerospace  flight  vehicles.  This  has  been  a  source 
of  serious  concern  since  there  are  no  analysis  techniques  available  that  have  predicted  LCO  in  an 
operational  aircraft.  There  have  been  some  semi-empirical  techniques  developed  to  correlate 
with  LCO  that  have  been  observed  in  flight,  and  these  are  useful  for  understanding  the  LCO  that 
has  occurred  (see  Ref  24).  However  these  techniques  are  not  as  satisfactory  for  the  design  of  a 
new  vehicle  or  the  substantial  modification  of  an  existing  one,  e.g.,  new  stores  to  be  carried  by 
an  aircraft.  In  this  regard,  it  should  be  noted  that  LCO  may  be  beneficial  as  well  as  detrimental. 
Without  the  nonlinearities  that  lead  to  LCO,  the  onset  of  flutter  may  lead  to  catastrophic  failure 
of  the  structure.  Hence  if  we  can  understand  and  predict  LCO,  perhaps  we  can  take  advantage  of 
these  nonlinearities  to  shape  more  favorable  responses  of  the  aircraft  leading  to  enhanced  safety 
and  performance. 

Sources  ofNonlinearities 

The  principal  sources  of  the  nonlinearities  essential  to  the  LCO  are  a  subject  of  current  debate 
among  the  experts  in  the  field.  The  candidate  sources  are  several: 

Fluid 

•  Shock  motions 

•  Separated  flow 

Structure 

•  Free-play 

•  Geometric,  e.g.,  a  nonlinear  relationship  between  strain  and  displacement 

•  Material,  e.g.,  dry  friction 

Also  there  is  a  further  distinction  between  a  static  versus  a  dynamic  nonlinearity.  An  important 
example  of  this  is  the  role  of  a  shock  wave  in  the  fluid.  If  a  shock  is  present,  then  its  creation  is 
the  result  of  a  dynamic  nonlinear  process.  However  once  a  steady  flow  is  established,  and  if  the 
airfoil  motion  is  sufficiently  small,  then  the  shock  motion  will  also  be  small  and  proportional  to 
the  airfoil  motion.  Hence  in  this  situation,  the  shock  itself  represents  a  nonlinear  static  (time 
independent)  equilibrium  and  the  motion  may  be  treated  as  a  dynamically  linear  perturbation 
about  the  mean  shock  position.  In  most  of  the  following  discussion,  we  assume  a  dynamically 
linear  model  of  the  shock  motion,  but  also  include  a  structural  (dynamic)  nonlinearity,  i.e., 
freeplay  in  the  connection  of  the  control  surface  to  the  airfoil. 

Airfoil  ysilk  Control  Surface  Free  play  Model 

A  sketch  of  the  configuration  is  shown  in  Fig  2.1.  It  is  a  conventional  typical  section  model 
except  that  the  spring  that  attaches  the  control  surface  to  the  airfoil  has  a  nonlinear  freeplay.  The 
elastic  restoring  torque  or  moment  provided  by  this  spring  is  shown  in  Fig  2.2  as  a  function  of 
control  surface  or  flap  rotation  angle,  J3.  The  freeplay  angle  is  8.  Note  that  when  f  is  less  than  8, 
there  is  zero  restoring  torque,  while  for  /?  greater  than  8,  the  spring  stiffness  is  the  nominal  value 
in  the  absence  of  freeplay.  The  freeplay  may  be  thought  of  as  creating  a  stiffness  or  (uncoupled) 
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natural  frequency  of  the  spring  that  varies  as  a  function  of  flap  amplitude.  This  interpretation  is 
shown  in  Fig  2.3.  Here  the  flap  uncoupled  natural  frequency  normalized  by  the  nominal  value  in 
the  absence  of  freeplay  is  shown  as  a  function  of  flap  amplitude,  /?,  normalized  by  the  freeplay 
angle,  S.  Note  that  given  a  certain  flap  amplitude,  there  is  a  corresponding  "equivalent"  flap 
frequency.  Of  course  for  a  linear  system  the  control  surface  or  flap  frequency  would  have  a 
fixed  value  independent  of  flap  amplitude. 


Figure  2.1  Airfoil  with  Control  Surface 


0.0 

Flap  Displacement,  p 

Figure  2.2  Restoring  Moment  As  a 
Function  of  Control  Surface  or  Flap 
Rotation  Angle 


Figure  2.3  Equivalent  Nonlinear  Flap  Frequency  Versus  Flap  Rotation 


Thus,  conceptually  and  computationally,  one  may  proceed  as  follows.  First,  one  determines  the 
neutrally  stable  motions  of  the  system  in  the  absence  of  freeplay  for  various  flap  frequencies 
from  zero  to  the  nominal  value.  Then  one  determines  the  corresponding  neutrally  stable 
nonlinear  limit  cycle  motion,  namely  the  flap  amplitude,  from  Fig  2.3. 
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This  approach  has  been  used  successfully  at  low  Mach  number  and  the  theoretical  results 
correlated  with  experiments  (see  Refs  25,  26). 

Computational  Fluid_  Dynamic^  (CFD)  Modeling  and_  its  Modal  Decomposition 

A  typical  CFD  model  is  very  large  in  terms  of  the  number  of  equations  required  to  be  solved. 
And  this  makes  such  models  problematical  for  aeroelastic  (and  some  other)  analyses.  For 
example,  the  CFD  model  used  in  the  present  work  is  based  upon  the  Euler  equations  of  fluid 
mechanics  and  has  a  spatial  grid  of  65  x  97  (6035)  mesh  points.  At  each  grid  point  there  are  four 
fluid  variables  to  be  determined.  Thus  the  CFD  model  per  se  has  about  25,000  unknowns  to  be 
determined  by  solving  25,000  equations.  This  is  a  doable  task  if  the  structural  motion  is  known. 
However  if  this  CFD  model  is  to  be  combined  with  a  set  of  structural  equations  of  motion,  and 
solutions  are  to  be  found  for  many  combination  of  structural  and  fluid  parameters,  then  the 
calculation  using  the  original  CFD  model  quickly  gets  out  of  hand.  Thus  the  search  for  an 
alternative  approach. 

Several  semi-empirical  methods  have  been  derived  to  address  the  computational  feasibility  issue. 
These  as  well  as  the  method  to  be  described  and  used  here  are  discussed  in  more  depth  in  Ref  1 1. 
Among  these  methods  are  variations  on  the  notion  of  a  Pade  approximant.  The  method  used 
here  is  based  upon  the  observation  that  virtually  all  CFD  models  can  be  thought  of  as  having  a 
modal  composition.  The  simplest  conceptual  set  of  modes  is  perhaps  the  fluid  or  aerodynamic 
eigenmodes  of  the  CFD  model,  and  these  modes  have  been  used  successfully  in  creating  reduced 
order  models  (ROM)  that  are  computationally  and  conceptually  attractive.  See  the  discussion  in 
Ref  1 1  and  the  forthcoming  Ref  27. 

However  it  turns  out  that  determining  the  aerodynamic  eigenmodes  of  a  large  CFD  model  is 
itself  a  challenging  task.  Hence  a  method  called  Proper  Orthogonal  Decomposition  (POD)  is 
employed  here  and  in  (Ref  1 1 ). 

For  the  present  analysis,  63  POD  modes  are  found  from  the  frequency  responses  (aerodynamic 
transfer  functions)  in  flap,  pitch  and  plunge  respectively  at  2 1  frequencies  at  each  Mach  studied 
using  the  original  CFD  model.  Based  upon  previous  experience,  one  might  use  an  even  smaller 
number  of  aerodynamic  modes  than  this.  However,  even  with  this  generous  number  of  modes, 
the  computations  described  below  were  all  done  in  a  few  days.  The  computational  grid  used  for 
the  CFD  model  is  shown  in  Fig  2.4. 
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(a)  Baseline  (b)  Grid  Motion 

Figure  2.4  NACA  0012  Grid:  65  x  97  Computational  Nodes,  Outer  Boundary  Radius=15c 


Linear  Instability  (Flutter) 


First,  consider  the  flutter  behavior  for  this 
system  in  the  absence  of  freeplay.  The 
stability  of  this  system  was  assessed  by 
constructing  a  root  locus  (migration  of  the  true 
aeroelastic  eigenvalues)  as  a  function  of  the 
nondimensional  airspeed  or  dynamic  pressure 
for  each  Mach  number.  The  usual  structural 
and  flow  parameters  are  defined  in  (Refs  24, 
25). 


A  representative  root  locus  result  is  shown  in 
Fig  2.5  for  M=0.80.  Root  locus  results  are 
available  for  all  Mach  numbers  used  to 
construct  the  flutter  boundary  which  is  shown 
in  Fig  2.6. 


Real  Frequency  Ratio,  Re(X/<Da) 

Figure  2.5  Aeroelastic  System  Eigenvalue 
Root-Loci:  M  =  0.80,  a  =  0.0  deg. 


Note  especially  that  the  "plunge"  aeroelastic  mode  has  a  root  that  for  low  "gain"  (or  flow 
velocity  or  dynamic  pressure)  moves  to  the  left  and  becomes  more  stable.  But  then  as  the  flow 
velocity  increases,  it  reverses  direction  and  moves  into  the  right  half  plane  becoming  unstable. 
And  then  at  even  higher  velocities,  it  moves  back  into  the  left  hand  plane  and  becomes  stable 
again.  However  by  then,  the  pitch  mode  has  moved  into  the  right  half  plane  and  become 


10 


unstable.  Hence  the  aeroelastic  system  remains  unstable  once  the  plunge  mode  becomes  again 
stable  at  this  velocity  and  Mach  number. 
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Figure  2.6  Mach  Number  Flutter  Trend:  a  =  0.0  deg. 


All  the  other  roots  in  this  figure  which  appear  to  originate  from  near  the  origin  are  essentially 
aerodynamic  roots,  and  these  roots  all  move  off  into  the  left  half  plane  indicating  they  are  always 
stable  and  increasingly  so  as  the  dynamic  pressure  increases.  As  the  Mach  number  becomes 
higher,  the  most  critical  root  may  change.  For  example,  at  M=0.85  the  pitch  root  becomes 
unstable  first,  and  for  M=0.90,  it  is  the  flap  root.  At  yet  higher  Mach  numbers,  no  roots  become 
unstable.  For  brevity  these  other  root  loci  are  omitted  here.  Taking  all  of  this  information  from 
the  root  loci  at  various-  Mach  numbers,  the  flutter  boundary  trend  with  Mach  number  can  be 
determined,  and  is  presented  in  Fig  2.6.a. 

There  are  several  interesting  features  to  this  flutter  boundary.  Up  to  M=0.80,  the  root-loci  are 
rather  similar,  and  it  is  always  the  plunge  root  that  is  critical  for  flutter.  Starting  at  M=0.80,  the 
pitch  root  also  shows  instability,  and  at  M=0.825  and  0.85,  it  is  most  critical  for  flutter.  At 
M=0.875  and  0.90,  the  flap  mode  is  most  critical  for  flutter,  and  for  M=0.925  to  at  least  M=l.l, 
no  flutter  is  observed  for  a  non-dimensional  flow  velocity  up  to  at  least  one.  The  corresponding 
frequencies  of  the  flutter  oscillations  are  shown  in  Fig  2.6.b.  Note  that  in  Fig  2.6.a  when  two 
data  points  are  shown  for  say  the  plunge  root  at  a  fixed  Mach  number,  the  lower  velocity  point  is 
when  flutter  begins,  and  the  higher  velocity  point  is  when  the  root  returns  to  the  stable  left  half 
plane  and  flutter  ceases  in  that  root.  Note  also  the  narrow  range  of  Mach  number  where  the 
change  in  flutter  mode  occurs.  Results  of  this  type  have  been  observed  in  experiments  where 
they  are  called  "chimneys"  (see  Ref  28). 

Limit  Cycle  Oscillations^ 

Now  the  freeplay  is  added  to  the  model  and  thus  LCO  may  occur.  As  is  perhaps  obvious  from 
physical  intuition,  when  freeplay  is  added,  the  stiffness  of  the  control  surface  freeplay  is  reduced 


11 


for  small  motions.  Hence  one  expects  limit  cycle  oscillations  to  occur  below  the  flutter 
boundary,  i.e.,  at  flow  velocities  less  than  those  shown  in  Fig  2.6. a.  Indeed  a  few  moments  of 
reflection  may  lead  one  to  expect  that  once  the  linear  flutter  boundary  shown  in  Fig  2.6.a  is 
exceeded,  then  exponentially  explosive  flutter  will  occur  when  the  nonlinearity  is  due  to 
ffeeplay.  That  is  found  to  be  the  case  as  shown  by  the  present  analysis  and  also  by  the  analysis 
and  experiments  of  Refs  25  and  26. 

The  LCO  results  shown  in  Fig  2.7  have  several  interesting  features.  First  of  all,  the  limit  cycle 
amplitude  is  normalized  by  the  ffeeplay  angle,  8.  The  theory  predicts  and  experiments  agree, 
(see  Ref  26),  that  when  the  results  are  normalized  in  this  manner,  they  are  universal.  That  is,  the 
limit  cycle  amplitude  is  proportional  to  the  ffeeplay  angle.  Secondly,  for  the  lowest  velocity  at 
which  LCO  may,  a  finite  disturbance  is  required  to  generate  LCO  at  this  lowest  velocity  and  for 
a  small  velocity  range  thereafter.  LCO's  for  any  disturbance,  no  matter  how  small,  will  only 
occur  when  the  flutter  velocity  for  a  flap  natural  ffequency  of  zero  is  exceeded,  at  a  somewhat 
higher  velocity,  about  0.22  in  Fig  2.7.a.  The  unstable  LCO's,  which  are  shown  along  with  the 
stable  LCO's  (those  that  are  observed  in  an  experiment),  provide  a  measure  of  the  level  of 
disturbance  required  to  initiate  the  LCO  at  the  lower  flow  velocities. 

The  results  Fig  2.7  are  typical  until  one  reaches  the  higher  transonic  Mach  numbers  where  linear 
theory  predicts  flutter  will  cease.  At  the  highest  Mach  number  considered  here  where  flutter  and 
LCO  may  occur,  M=0.90,  the  LCO  has  a  somewhat  different  character.  Again  the  LCO  is  first 
encountered  at  the  minimum  velocity  at  which  flutter  will  occur  over  the  range  of  flap 
frequencies.  But  now  the  corresponding  flap  frequency  is  zero.  Moreover,  when  the  flow 
velocity  increases  to  higher  values,  there  are  two  stable  limit  cycles.  The  nature  of  the 
disturbances  to  the  system  would  determine  which  of  these  two  LCO  would  be  observed  in  a 
wind  tunnel  experiment  or  in  flight.  For  brevity,  the  M=0.9  results  are  not  shown. 
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Flutter  Velocity,  VF-Uf/p  co0b 


(a)  LCO  Flap  Rotation 


(b)  LCO  Frequency  Ratio 


Figure  2.7  Limit  Cycle  Behavior:  M  =  0.80,  a  =  0.0  deg. 

For  additional  details  of  the  work  presented  in  this  Chapter,  please  see  the  paper  (Ref  14) 
included  in  Appendix  A. 
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CHAPTER  3 


NONLINEAR  INVISCID  AERODYNAMIC  EFFECTS  ON  TRANSONIC 
DIVERGENCE,  FLUTTER  AND  LIMIT  CYCLE  OSCILLATIONS 


Summary 


Aerodynamic  nonlinearities  may  be  give  rise  to  LCO,  and  these  may  be  either  stable  (favorable) 
or  unstable(unfavorable).  An  example  of  the  former  is  shown  here  as  the  nonlinear  counterpart 
of  classical  linear  aeroelastic  divergence.  An  example  of  the  latter  is  also  shown  here  as  the 
nonlinear  counterpart  of  single-degree-of-freedom  pitch  flutter.  Future  work  will  be  directed 
toward  the  study  of  the  nonlinear  counterpart  of  classical  bending/torsion  flutter  where  similar 
methods  may  be  used. 


Introduction 


Limit  cycle  oscillations  (LCO)  in  aeroelastic  systems  appear  to  be  more  prevalent  in  transonic 
flow  than  in  subsonic  flow.  Hence  it  has  been  thought  that  at  least  for  some  configurations  the 
source  of  the  nonlinearity  that  leads  to  LCO  is  in  the  aerodynamic  flow.  See  Ref  29  and  26  for  a 
discussion  of  structural  nonlinearities  relevant  to  LCO.  Thus,  we  consider  the  effects  of 
nonlinearities  arising  from  inviscid  transonic  aerodynamics.  The  principal  physical  effect  of 
interest  is  the  relatively  large  motion  of  the  shock  wave  as  the  amplitude  of  say  the  pitch  motion 
of  the  airfoil  becomes  sufficiently  large.  This  in  turn  leads  to  a  movement  of  the  center  of 
pressure  with  amplitude.  Hence  one  expects  to  see  an  effect  of  amplitude  on  the  neutrally  stable 
motions  that  may  occur.  Moreover  this  may  lead  to  limit  cycle  motions  rather  than  the 
catastrophic  exponentially  growing  oscillations  predicted  by  time  linearized  models.  The  latter 
models  capture  the  effect  of  the  mean  position  of  the  shock  and  small  shock  motions  about  this 
mean  position  by  assuming  the  shock  motion  is  dynamically  linear,  i.e.,  the  shock  motion  is 
proportional  to  the  airfoil  motion.  This  is  not  true  for  dynamically  nonlinear  aerodynamic 
models  that  allow  for  larger  and  more  general  shock  motions  including  the  possible  appearance 
and  disappearance  of  a  shock  during  a  cycle  of  airfoil  motion.  The  latter  is  our  concern  here. 

We  will  consider  two  distinct  aeroelastic  phenomena,  divergence  and  flutter,  and  their  associated 
limit  cycle  oscillations.  To  keep  the  discussion  focussed  on  the  fundamental  physical 
phenomena,  and  to  ease  the  interpretation  of  the  inherently  complex  phenomena,  only  a  single 
structural  degree  of  freedom  will  be  studied.  However  the  aerodynamic  model  will  be  a  state-of- 
the-art  computational  fluid  dynamics  (CFD)  based  upon  the  Euler  equations  of  nonlinear, 
rotational  inviscid  aerodynamic  theory. 

Here  we  emphasize  that  the  solution  technique  is  for  a  large  system  of  ordinary  differential 
equations  in  time,  which  represents  the  time  variation  of  the  fluid  unknowns  at  each  spatial  grid 
point  in  the  CFD  model.  The  unknowns  are  four  in  number  at  each  grid  point  for  a  two- 
dimensional  Euler  flow  and,  for  example,  could  be  density,  the  two  scalar  components  of 
momentum,  and  the  total  energy  at  each  grid  point.  The  present  CFD  model  has  about  17,000 
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total  flow  variable  unknowns,  and  therefore  an  efficient  solution  method  is  imperative  to  carry 
out  the  studies  reported  here. 

Harmonic^  Balance^  Solution  in  the  Frequency  Domain 

The  pioneering  work  of  Ueda  and  Dowell,  (Ref  30)  and  Lan  and  his  colleagues,  (Ref  31)  should 
be  recalled.  Ueda  and  Dowell  used  a  describing  function  technique  whereby  the  dominant 
harmonic  was  extracted  from  a  time  marching  CFD  model,  LTRAN2,  using  both  indicial  and 
harmonic  motions  of  the  airfoil.  They  considered  a  two  degree  of  freedom  typical  airfoil  section. 
Lan  et  al  used  the  method  of  harmonic  balance  to  study  the  unsteady  transonic  aerodynamics  for 
flutter  and  limit  cycle  oscillation  prediction.  In  their  work,  they  used  the  transonic  small 
disturbance  potential  flow  model,  as  did  Ueda  and  Dowell,  and  only  considered  a  single 
harmonic.  In  the  present  work,  we  employ  the  Euler  equations  of  fluid  dynamics  and  also  retain 
multiple  harmonics  in  the  aerodynamic  model.  It  is  found  that  using  several  harmonics  improves 
the  theoretical  prediction  of  the  aerodynamic  forces.  However  in  the  aeroelastic  analysis,  when 
the  fluid  and  structural  models  are  coupled,  only  a  single  harmonic  is  used.  The  effects  of  higher 
harmonics  on  this  single  harmonic  are  retained  as  they  are  found  to  be  significant  in  the  fluid 
model. 

The  Aeroelastic^  System  and_  Its  Solution 

The  structural  equation  of  motion  is  a  simple  single  degree  of  freedom  model  in  pitch.  See  Fig 
3.1  for  a  depiction  of  the  airfoil  and  the  CFD  grid  used  in  the  numerical  calculations.  By 
carefully  selecting  the  pitch  axis  and  mass  ratio,  we  can  insure  that  the  system  will  either 
undergo  classical  linear  aeroelastic  divergence  or  flutter.  Divergence  can  occur  when  the 
aerodynamic  "negative"  stiffness  overcomes  the  structural  stiffness,  while  flutter  may  occur 
when  the  aerodynamic  negative  damping  overcomes  the  structural  damping.  As  will  be  shown, 
each  of  these  classical  linear  aeroelastic  phenomena  has  a  distinctively  different  limit  cycle  or 
nonlinear  behavior. 


Figure  3.1  NACA  64A010A  Computational  Grid 
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The  Mach  number  for  these  studies  is  M=0.80 
and  a  NACA  64A010A  airfoil  is  considered. 
The  NACA  64A010A  is  a  symmetric  (10.6% 
thickness  ratio)  variant  of  the  "Ames"  AGARD 
156  benchmark  section.  The  elastic  axis  is 
considered  at  the  mid-chord.  Employed  for  the 
CFD  calculations  is  an  "0"-type 
computational  mesh  with  65  x  65  radial  and 
circumferential  nodes  that  has  an  outer 
boundary  radius  of  10  chord  lengths.  The 
center  of  pressure  (xcp)  as  a  function  of  static 
angle  of  attack  is  shown  Fig  3.2  where  it  is 
seen  the  center  of  pressure  moves  from  32% 
chord  to  40%  chord  as  the  angle  of  attack 
varies  from  0.0  to  5.0  degrees.  This  is  a  key 
characteristic  of  the  flow  field  for  LCO. 


Airfoil  Center  of  Pressure,  x_/c 

cp 

Figure  3.2  Center  of  Pressure  Variation 
with  Angle  of  Attack:  NACA  64A010A 
Airfoil  Section,  M  =  0.80 


Lineai '  and  Nonlinear  Divergence 


Divergence  is  perhaps  the  simpler  of  the  two  phenomena  since  by  definition  it  is  time 
independent,  i.e.,  we  are  dealing  with  a  static  linear  instability  and  its  nonlinear  counterpart.  In 
this  case,  the  structural  equation  of  motion  becomes  an  equation  of  static  equilibrium  And  for 
the  aerodynamic  model,  we  only  need  to  determine  the  lift  and  moment  about  some  appropriate 
axis  as  a  function  of  angle  of  attack.  For  small  angle  of  attack,  we  will  recover  the  classical 
linear  aeroelastic  divergence  phenomena.  But  the  question  is,  what  are  the  effects  of  the 
nonlinearity? 


Qualitatively  one  can  anticipate  the  effect  of  the  aerodynamic  nonlinearity  by  examining  the 
aerodynamic  moment  variation  with  angle  of  attack.  A  necessary  condition  for  divergence  to 
occur  is  that  the  aerodynamic  moment  be  positive  in  the  same  direction  as  the  twist  angle 
Moreover,  if  the  nonlinear  aerodynamic  model  predicts  a  moment  less  in  magnitude  than  that 
predicted  by  linear  aerodynamic  theory,  the  effect  of  the  nonlinearity  will  be  to  stabilize  the 
divergence.  And  vice  versa  if  the  nonlinear  theory  predicts  an  increase  in  aerodynamic  moment 
over  that  given  by  linear  theoiy.  Hence  by  examining  the  slope  of  the  moment  vs.  angle  of 

attack  curve  with  increasing  angle  of  attack,  we  will  know  whether  the  effect  of  the  nonlinearity 
is  favorable  or  unfavorable.  J 


n  t  e  example  below,  the  effect  is  favorable.  That  is,  once  the  divergence  dynamic  pressure  for 
a  small  angle  of  attack  is  exceeded  (this  is  the  classical  linear  aeroelastic  divergence  dynamic 
pressure),  then  the  angle  of  twist  of  the  pitch  spring  remains  finite  and  smoothly  increases  from 
zero  beyond  the  divergence  dynamic  pressure. 
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See  Fig  3.3  where  the  angle  of  twist  is  plotted 
vs.  the  non-dimensional  dynamic  pressure. 
Also  shown  are  results  with  an  initial  angle  of 
attack.  In  this  latter  case,  there  is  some  twist 
over  the  full  range  of  dynamic  pressure. 
Indeed  even  if  the  initial  angle  of  attack  is  only 
a  few  degrees,  it  would  be  difficult  to  detect 
the  classical  divergence  dynamic  pressure 
experimentally  for  this  example.  For  readers 
who  have  studied  buckling  of  systems  in  the 
presence  of  imperfections  (e.g.  beams,  plates 
or  shells  with  initial  curvature),  this  behavior 
will  be  familiar. 


Nondimensional  Dynamic  Pressure.  X=q_c2/Ka  (radians) 
Figure  3.3  Divergence  and  Post-Divergence 
of  an  Airfoil  Including  Transonic  Nonlinear 
Inviscid  Aerodynamics:  NACA  64A010A 
Airfoil  Section,  M  =  0.80,  Elastic  Axis 
Location,  a  =  e/b  =  0.0 


In  this  example,  recall  the  center  of  pressure  moves  from  32%  chord  at  low  angles  of  attack  to 
40%  chord  at  5.0  degrees  angle  of  attack.  This  is  the  principal  reason  for  the  stabilizing  effect  of 
nonlinear  aerodynamics  on  the  post-divergence  condition.  Had  the  change  of  the  slope  of  the 
aerodynamic  moment  curve  been  in  the  opposite  direction,  then  the  angle  of  twist  vs.  dynamic 
pressure  curve  would  have  bent  the  other  way.  That  is,  for  dynamic  pressures  below  the  classical 
divergence  dynamic  pressure,  there  would  be  non-trivial  (non-zero)  twist  angles  that  represent 
possible  static  nonlinear  equilibrium  solutions.  Intuitively  one  recognizes  that  these  latter 
solutions  would  themselves  be  unstable,  i.e.,  such  results  would  be  interpreted  physically  as  the 
magnitude  of  the  disturbance  required  to  generate  non-trivial  twist  at  dynamic  pressures  below 
the  classical  divergence  dynamic  pressure.  In  our  studies  to  date,  only  the  stable  nonlinear  effect 
has  been  observed  for  statically  divergent  systems.  However,  this  is  not  to  say  that  unstable 
nonlinear  divergence  systems  may  not  be  encountered  for  some  other  parameter  combinations. 

Of  course,  divergence  is  a  very  special  case  of  nonlinear  aeroelasticity  as  it  is  for  linear 
aeroelasticity,  because  the  frequency  of  oscillation  is  zero  when  divergence  and  post-divergence 
occurs.  Thus  we  now  turn  to  an  oscillatory  case. 


Flutter  and  Associated  LCO 


Now  consider  single-degree-of-ffeedom  flutter  in  pitch.  Here  the  classical  flutter  arises  from  a 
negative  damping  in  the  aerodynamic  moment  beyond  a  certain  reduced  frequency.  However  the 
reduced  frequency  at  which  the  aerodynamic  damping  moment  becomes  negative  increases  as 
the  angle  of  pitch  oscillation  increases.  Hence  the  reduced  velocity  decreases  as  the  angle  of 
pitch  increases,  which  suggests  that  this  will  lead  to  an  unstable  LCO  as  indeed  it  does. 

In  the  example  considered,  we  have  moved  the  elastic  axis  to  20%  chord  to  preclude  divergence 
and  to  induce  flutter. 
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It  should  be  emphasized  that  in  the  present  analysis,  we  are  using  a  single  harmonic  to  represent 
the  pitch  oscillation.  However  in  the  calculation  of  the  aerodynamic  moment,  we  include  up  to 
three  harmonics  to  determine  the  effect  of  higher  harmonics  on  the  first  harmonic  of  the 
aerodynamic  moment.  It  turns  out  that  the  effect  of  the  third  harmonic  is  negligible.  Indeed,  if 
one  only  retains  a  single  harmonic  in  the  aerodynamic  analysis,  the  results  are  qualitatively 
correct  and  have  fair  quantitative  accuracy. 

Results  for  the  first  harmonic  for  the  moment  about  the  pitch  or  elastic  axis  are  shown  in  Fig  3.4. 
These  results  are  for  two  harmonics  retained  in  the  aerodynamic  analysis.  Note  that  the  results  at 
a  reduced  frequency  of  zero  were  those  used  in  the  divergence  analysis  discussed  previously.  Of 
course,  a  transformation  of  the  pitch  axis  is  used  for  the  divergence  analysis. 


(a)  Real  Unsteady  Moment  (b)  Imaginary  Unsteady  Moment 

Figure  3.4  Unsteady-Lift  and  Moment  for  Various  Pitch  Amplitudes:  NACA  64A010A 
Airfoil  Section,  M  =  0.8,  a0  =  0.0  deg.,  Elastic  Axis  Location,  a  =  e/b  =  -0.6,  Two  Harmonics 
Employed  in  Harmonic  Balance  Expansion 

While  structural  damping  is  readily  included  in  the  analysis,  it  will  be  helpful  to  understand  the 
essence  of  the  results  by  first  considering  the  solution  for  zero  structural  damping. 

Zero  Structural  Damping 

In  this  case,  a  neutrally  stable  oscillation  will  occur  when  the  imaginary  part  of  the  aerodynamic 
moment  becomes  zero.  This  will  occur  at  some  reduced  frequency  for  a  particular  angle  of  pitch 
oscillation  (and  other  parameters  fixed  such  as  Mach  number).  Then  one  can  solve  for  the 
frequency  of  this  neutrally  stable  oscillation.  For  sufficiently  small  motions,  this  is  the  flutter 
solution;  for  larger  motions,  we  determine  a  limit  cycle  oscillation.  The  solution  procedure  then 
is  to  select  an  amplitude  of  oscillation,  determine  the  reduced  frequency  at  which  the  imaginary 
part  of  the  aerodynamic  moment  is  zero  from  Fig  3.4,  and  then  determine  the  frequency  of  the 
oscillation.  Note  this  is  essentially  the  same  computational  procedure  as  for  a  classical  flutter 
solution,  except  now  the  reduced  frequency,  the  frequency  of  oscillation,  and  the  reduced 
velocity  are  all  functions  of  the  pitch  amplitude.  It  should  be  noted  however  that  just  because  the 
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imaginary  part  of  the  aerodynamic  moment  vanishes  (i.e.,  the  aerodynamic  damping  becomes 
zero),  that  alone  does  not  insure  that  a  neutrally  stable  oscillation  will  occur.  This  is  because  the 
frequency  determined  must  be  physically  possible. 


Large  Pitch  Moment  of  Inertia  Now  if  the 
mass  ratio  or  moment  of  inertia  is  larger,  a  not 
uncommon  circumstance,  then  the  flutter  or 
LCO  frequency  is  simply  equal  to  the 
structural  pitch  natural  frequency.  With  this 
approximation,  the  results  of  Fig  3.5  are 
obtained  for  both  zero  and  non-zero  structural 
damping.  Note  that  the  curves  bend  to  the  left 
which  is  indicative  of  an  unstable  LCO.  That 
is,  these  results  are  to  be  interpreted  as  the 
amplitude  of  a  disturbance  required  to  initiate 
explosive  flutter  below  the  classical  flutter 
velocity  for  this  single-degree-of-freedom  pitch 
oscillation. 


Effects  of  Finite  Pitch  Moment  ofjnertia 

For  general  values  of  moment  of  inertia  and  structural  damping,  the  solution  algorithm  proceeds 
as  follows.  First  select  a  Mach  number  and  pitch  axis,  and  for  a  range  of  pitch  amplitudes, 
determine  the  first  harmonic  of  the  aerodynamic  moment  (including  higher  harmonics  of  the 
aerodynamic  model  and  their  effect  on  the  fundamental  harmonic).  Then  for  a  given  pitch 
amplitude,  choose  a  reduced  frequency  and  determine  the  flutter  or  LCO  oscillation  frequency. 
This  frequency  will  be  proportional  to  the  pitch  structural  frequency,  of  course.  With  the  flutter 
or  LCO  frequency  determined,  and  the  reduced  frequency  selected,  one  then  knows  the  flow 
velocity  corresponding  to  the  chosen  pitch  amplitude.  Finally,  determine  the  structural  damping 
value  necessary  to  give  a  neutrally  stable  flutter  or  limit  cycle  oscillation.  From  this  perspective, 
the  flutter  condition  is  simply  the  neutrally  stable  motion  that  may  exist  at  small  angles  of  twist, 
and  the  LCO  are  the  neutrally  stable  oscillations  that  may  exist  when  the  pitch  amplitude  is 
finite.  Of  course  the  flutter  or  LCO  may  become  unstable  when  it  is  perturbed  (e.g.,  by 
perturbations  in  the  amplitude  of  oscillation),  and  this  is  indeed  the  case  in  the  example  treated 
here. 

Up  to  this  point,  we  have  assumed  that  the  pitch  moment  of  inertia  is  well  above  its  asymptotic 
value.  Hence  the  flutter  frequency  is  the  same  as  the  structural  natural  pitch  frequency. 
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Reduced  Velocity,  U/<oac 
Figure  3.5  LCO  Amplitude  Versus 
Reduced  Velocity:  NACA  64A010A  Airfoil 
Section,  M  =  0.80,  a  =  0.0  deg.,  Elastic  Axis 
Location,  a  =  e/b  =  -0.6. 


18 


Now  we  consider  the  more  general  case  and  a 
range  of  pitch  inertias  such  that  the  flutter 
frequency  is  no  longer  precisely  equal  to  the 
structural  natural  frequency  in  pitch.  Results 
are  shown  for  non-dimensional  pitch  inertias  of 
200,100,50,37.5  and  25  in  Fig  3.6.  These  are 
for  LCO  amplitude  versus  reduced  velocity. 
The  asymptotic  or  large  pitch  inertia  results  are 
also  shown  for  reference. 

As  expected,  for  sufficiently  large  pitch  inertia, 
say  greater  than  200,  the  asymptotic  results  are 
good  approximations.  However  for  pitch 
inertias  less  than  1 00,  the  results  show  a  more 
sensitive  dependence  on  pitch  moment  of 


Reduced  Velocity,  U/a>0c 
Figure  3.6  LCO  Amplitude  Versus 
Reduced  Velocity  for  Various  Pitch  Inertias: 
NACA  64A010A  Airfoil  Section,  M  =  0.80,  a 
=  0.0  deg.,  Elastic  Axis  Location,  a  =  e/b  =  - 


inertia.  For  sufficiently  small  pitch  moment  of 
inertia,  of  course,  no  flutter  or  LCO  is 
possible. 


For  additional  details  of  the  work  presented  in  this  Chapter,  please  see  the  paper  (Ref  15) 
included  in  Appendix  B. 
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CHAPTER  4 


THREE  DIMENSIONAL  TRANSONIC  AEROELASTICITY  USING 
PROPER  ORTHOGONAL  DECOMPOSITION  BASED  REDUCED  ORDER 

MODELS 


Summary 

The  POD/ROM  method  has  been  demonstrated  for  the  flutter  analysis  of  a  three-dimensional 
transonic  wing  configuration.  We  have  shown  that  the  number  of  ROM  DOF's  necessary  to 
create  accurate  models  is  on  the  order  of  a  few  dozen  as  is  the  case  in  two-dimensions.  We  have 
also  shown  that  it  is  unnecessary  to  compute  a  completely  new  ensemble  of  solution  snapshots 
based  on  the  vibratory’  mode  shapes  for  each  new  structural  configuration  that  might  be  under 
consideration.  One  can  simply  compute  a  set  of  snapshots  based  on  some  basic  wing  motions  at 
a  number  of  frequencies.  Then  snapshots  only  at  the  end  points  of  the  frequency  range  of 
interest  need  to  be  computed  for  the  specific  mode  shapes  of  the  configuration  of  interest.  These 
end  point  snapshots  ‘‘lock  in  ”  the  unsteady  fluid  dynamic  characteristics  for  the  particular  mode 
shapes,  and  the  simple  motion  snapshots  then  act  to  resolve  the  dominant  dynamics  of  the  flow 
throughout  the  full  frequency  range  of  interest. 


Introduction 

Here,  we  demonstrate  how  the  recently  devised  proper  orthogonal  decomposition  (POD)  based 
reduced  order  modeling  (ROM)  technique  (Refs  11,  32)  can  be  used  to  model  unsteady 
aerodynamic  and  aeroelastic  characteristics  of  three-dimensional  transonic  wing  configurations. 
Although  transonic  Euler  flows  are  considered  in  Refs  1 1  and  32,  the  initial  demonstrations  of 
the  POD/ROM  method  as  presented  in  these  references  are  for  two-dimensional  flow  and  two 
structural  degree-of-ffeedom  airfoil  configurations.  Also  in  Ref  33,  an  application  of  the 
POD/ROM  technique  to  the  well  known  vortex  lattice  method  has  been  presented. 

In  extending  the  POD/ROM  technique  to  three-dimensions,  two  primary  issues  have  been  of 
concern.  First,  the  size  of  the  computational  fluid  dynamic  (CFD)  model  will  in  general  be  at 
least  an  order  of  magnitude  greater  than  for  two-dimensions.  Whereas  a  typical  CFD  model  for  a 
realistic  two-dimensional  configuration  might  have  on  the  order  of  10's  or  even  100's  of 
thousands  of  degrees  of  freedom  (DOF),  a  CFD  model  for  a  three-dimensional  configuration 
might  easily  have  on  the  order  of  at  least  hundreds  of  thousands  if  not  millions  or  more  DOF's. 
In  two-dimensions,  we  have  found  that  very  accurate  ROM's  with  on  the  order  of  only  a  few 
dozen  DOF's  can  be  devised  using  the  POD/ROM  methodology.  A  first  issue  to  address  has  thus 
been  whether  or  not  in  three-dimensions  one  can  also  generate  accurate  ROM's,  which  require  at 
most  a  few  dozen  DOF's. 

The  second  concern  is,  for  any  variation  of  the  structural  properties  of  a  given  wing  under 
consideration,  will  a  completely  new  ensemble  of  solution  vector  “snapshots”  have  to  be 
computed  in  order  to  devise  an  accurate  POD/ROM.  A  basic  aspect  of  the  POD/ROM  method  is 
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that  an  ensemble  of  solution  vectors  is  first  assembled  by  computing  unsteady  CFD  solutions  at  a 
number  of  discrete  frequencies  within  a  frequency  range  of  interest  for  the  unsteady  structural 
motions  that  are  also  of  interest.  In  two  dimensions,  this  step  is  relatively  straight  forward  since 
one  only  has  to  consider  a  few  possible  motions,  e.g.,  pitch  and  plunge. 

In  three-dimensions  however,  the  wing  vibratory  mode  shapes  will  be  different  for  each  different 
structural  configuration  of  a  given  wing.  There  can  in  fact  be  an  infinite  number  of  unsteady 
motions  (or  at  least  a  substantial  number  of  motions  equivalent  to  the  number  of  DOF's  of  the 
discrete  structural  model).  Thus  the  second  concern  about  extending  the  POD/ROM  to  three- 
dimensions  has  been  whether  or  not  it  is  necessaiy  to  compute  a  completely  different  ensemble 
of  solution  snapshots  for  every  possible  structural  configuration.  For  example,  say  one  computes 
solution  snapshots  for  a  given  wing  configuration  based  on  the  wing's  particular  vibratory  modes 
shapes  in  order  to  develop  a  POD/ROM  to  model  the  configuration's  aeroelastic  characteristics. 
Then  the  question  is,  if  the  structural  make-up  of  the  wing  changes,  does  one  have  to  compute  a 
whole  new  ensemble  of  solution  snapshots  for  the  same  wing,  but  for  the  different  set  of 
vibratory  modes  shapes. 

Fortunately  in  addressing  these  two  issues,  we  have  found  that  accurate  POD/ROM’s  with  just  a 
few  dozen  degrees  of  freedom  can  in  fact  be  created  for  a  realistic  transonic  three-dimensional 
configurations.  This  is  true  even  though  in  the  model  problem  to  be  shown  subsequently,  the 
CFD  model  is  easily  an  order  of  magnitude  larger  than  anything  we  have  previously  studied  in 
two-dimensions.  Furthermore,  we  have  discovered  that  a  “fundamental”  ensemble  of  solution 
snapshots,  based  on  wing  motions  that  need  not  be  related  to  the  structural  modes  under 
consideration,  can  be  assembled  as  a  first  step.  Accurate  POD/ROM's  for  a  given  wing 
configuration  can  then  be  created  by  simply  adding  to  this  “fundamental”  ensemble,  the 
snapshots  corresponding  to  actual  wing  structural  modal  motions  solely  at  the  frequencies 
corresponding  to  the  end  points  of  the  frequency  range  of  interest.  In  general,  these  two 
snapshots  prove  to  be  sufficient  to  “lock  in”  the  conditions  corresponding  to  the  particular 
structural  motion,  and  indeed  the  fundamental  ensemble  of  solution  snapshots  is  sufficient  to 
reveal  the  unsteady  dynamics  of  the  fluid  dynamic  model.  The  fundamental  ensemble  of 
snapshots  can  be  used  again  and  again  even  as  the  structural  mode  change,  and  thus  the 
computational  cost  of  having  to  compute  an  entirely  new  snapshot  ensemble  for  every  new 
structural  configuration  is  greatly  reduced. 

POD/ROM_  Methodology 


In  the  following,  we  will  be  considering  inviscid  three-dimensional  Euler  flows.  More 
specifically,  linearized  (about  some  nonlinear  background  steady  flow)  unsteady  frequency- 
domain  CFD  solutions  to  the  Euler  equations  are  computed.  The  POD/ROM  procedure  can  be 
considered  as  a  “wrapper”  around  any  typical  CFD  method,  and  the  CFD  method  we  have 
employed  for  the  present  analysis  is  a  variant  of  Ni's  (Ref  34),  approach  to  the  standard  Lax- 
Wendroff  method.  The  frequency  domain  CFD  method  in  effect  represents  a  linear  system 
formulation  of  the  unsteady  fluid  dynamic  model,  i.e. 


Aq  =  -B£ 


(4.1) 
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where  q  is  an  N  dimensional  vector  (N  is  the  number  of  mesh  points  time  the  number  of 
dependent  flow  variables)  of  the  unknown  flow  variables  a  each  mesh  point  in  the  CFD  domain, 
and  £"is  the  L  dimensional  vector  (L  is  the  number  structural  mode  shapes)  of  modal  coordinates 
for  the  structural  model.  A  is  the  Nx  N  fluid  dynamic  influence  matrix,  and  B  is  the  Nx  L  matrix 
which  relates  the  flow  solver  boundary  conditions  to  each  particular  mode  shape.  Both  A  and  B 
are  functions  of  the  background  flow  and  unsteady  frequency  co.  The  structural  equations  for  the 
wing  configuration  being  modeled  within  the  flow  can  be  written  as 


DC  =  ~Cq  (4.2) 

where  D  represents  the  L  x  L  structural  influence  matrix  (i.e.,  D  =  -co2 M  +  K  where  M  and  K 
are  the  generalized  mass  and  stiffness  matrices,  and  C  is  the  L  x  M  matrix  which  represents  the 
discrete  integration  used  to  obtain  the  generalized  forces  associated  with  each  modes  shape  based 
on  the  unsteady  flow  q.  When  Eqs  (4. 1 )  and  (4.2)  are  put  together, 


A 

C 


0 

0 


(4.3) 


The  resulting  Eq  (4.3)  is  a  fully  coupled  aeroelastic  system  of  equations,  which  for  nontrivial  q 
and  C  represents  an  eigenvalue  problem  with  co  being  the  eigenvalue.  Any  eigenvalues  with  a 
positive  real  part  imply  the  aeroelastic  system  is  unstable. 

The  problem  with  constructing  and  solving  this  eigenvalue  problem  is  that  A  is  simply  too  large 
for  realistic  configurations.  As  mentioned  in  the  introduction,  N  can  easily  be  on  the  order  of 
10,000  or  100,000  for  two-dimensional  configurations,  and  on  order  of  100,000  to  1,000,000  or 
even  more  for  three  dimensional  configurations.  For  such  large  cases,  even  attempting  to  set  up 
A  is  well  beyond  the  memory  limits  of  todays  largest  computers. 

The  basic  premise  of  the  POD/ROM  methodology  is  that  we  assume  the  unknown  flowfield 
solution  vector  q  can  be  expressed  as  a  Ritz  type  expansion  of  the  form 


q*E**.M 


K«N 


(4.4) 


where  Ck  *s  a  generalized  coordinate  sometimes  referred  to  as  an  augmented  aerodynamic  state 
variable,  and  <pk  is  the  corresponding  Ritz  vector.  Eq  (4.4)  can  also  be  written  in  matrix  form  as 


q  =  <t>£ ,  where  O  = 

-  -C  ' 

. -  ^  ' 

1 ' 
..  <t>K 

and  C  - ' 

£ 

r 

LI  1 

1 . 

£k. 
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Here,  ®  is  an  N  x  K  matrix  whose  k’th  column  is  the  shape  vector  <f>k ,  and  £  is  the  K 
dimensional  vector  of  augmented  aerodynamic  state  variables  . 

A  reduced-order  representation  of  the  fluid  dynamic  and  aeroelastic  systems  can  be  formulated 
by  substituting  Eq  (4.5)  into  Eq  (4.1)  and/or  (4.2)  and  pre-multiplying  by  the  Hermitian 

transpose  (® H )  of  ® ,  i.e. 


(D',AO^(PwB(  or  = 
and 


{® 


H 


dJU  j 
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(4.6) 


(4.7) 


If  the  Ritz  approximation  is  a  good  one,  (K  «  N),  and  Eqs  (4.6)  and  (4.7)  will  represent  much 
smaller  systems  that  can  readily  be  solved  using  conventional  eigenvalue  techniques. 


The  next  question  becomes  what  are  good  choices  for  the  Ritz  vectors  (f>k  that  will  in  fact  result 
in  good  Ritz  approximations.  Previous  studies  as  detailed  in  Refs  1 1  and  32  have  demonstrated 
that  shape  vectors  derived  via  the  proper  orthogonal  decomposition  technique  (see  for  instance 
Refs  35,  36,  and  37)  are  an  excellent  source.  For  the  sake  of  brevity,  the  of  details  are  omitted 
here,  but  a  discussion  of  how  the  shapes  are  derived  can  be  found  in  Refs  1 1  and  32.  The  basic 
premise  behind  their  formulation  is  that  a  number  solution  "snapshots"  are  directly  computed  for 
a  number  of  discrete  frequencies  and  unsteady  structural  motions  of  interest.  From  this  ensemble 
of  solution  vectors,  the  POD  shapes  are  easily  derived  by  solving  a  small  (the  size  of  the  number 
of  snapshots)  eigenvalue  problem.  The  first  few  POD  modes  describe  the  most  dominant 
dynamic  characteristics  of  the  fluid  dynamic  system,  and  as  such,  the  POD  shapes  have  proven 
to  be  an  excellent  set  of  Ritz  vectors  for  fluid  dynamic  and/or  aeroelastic  models. 


Model  Problem 


The  configuration  under  consideration  is  the  AGARD  model  445.6  wing  (Refs  38,  39).  This  is  a 
45  degree  quarter  chord  swept  wing  using  the  NACA  64A004  airfoil  section  that  has  an  aspect 
ratio  of  3.3  (for  the  full  span)  and  a  taper  ratio  of  2/3.  Fig  4.1  illustrates  the  computational  mesh 
employed  for  this  configuration.  The  grid  is  an  “O-O”  topology  that  employs  49  computational 
nodes  about  each  airfoil  section,  33  nodes  normal  to  the  wing,  and  33  nodes  along  the  semispan. 
The  outer  boundary  of  the  grid  extends  five  semispans  from  the  midchord  of  the  wing  root 
section.  The  particular  structural  configuration  of  the  wing  is  referred  to  as  the  2.5  ft.  weakened 
model  3  (again  see  Refs  38  and  39). 
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(a)  Wing  surface  and  symmetry  plane  grids 


(b)  Outer  boundary  grid 


Figure  4.1  AGARD  445.6  Wing  Grid  Topology 

We  have  examined  the  computed  wing  surface  and  symmetry  plane  steady  flow  pressure 
contours  for  ffeestream  Mach  numbers  of  0.960  and  1.141.  A  relatively  weak  shock  can  be  seen 
at  the  trailing  edge  for  M=0.960.  This  shock  appears  to  get  stronger  at  M=1.141.  The  wing 
section  is  quite  thin  4%,  so  a  strong  shock  is  not  really  expected.  Comparing  contours,  our 
flowfields  look  very  comparable  to  those  of  Lee-Rausch  and  Batina  (Ref  40),  although  they 
employed  a  much  larger  mesh.  We  have  examined  the  computed  wing  surface  and  symmetry 
plane  steady  flow  pressure  contours  for  the  Mach  numbers  of  0.960  and  1.141.  A  relatively 
weak  shock  can  be  seen  at  the  trailing  edge  for  M=0.960.  This  shock  appears  to  get  stronger  at 
M=1 . 141 .  The  wing  section  is  quite  thin  (4%),  so  a  strong  shock  is  not  really  expected. 
Comparing  contours,  our  flowfields  look  very  comparable  to  those  of  Lee-Rausch  and  Batina 
(Ref  40),  although  they  employed  a  much  larger  mesh. 

Computational  Steps  for  Flutter  Analysis 

The  first  step  is  to  compute  the  “snapshots”  for  one  Mach  number  and  several  structural  modes, 
say  five.  These  are  overnight  runs,  so  the  wall  clock  time  is  about  24  hours  to  do  this.  Moreover 
we  have  shown  that  if  the  structural  modes  change,  then  the  computation  of  all  of  the 
“snapshots”  does  not  need  to  be  repeated.  This  reduces  the  number  of  new  “snapshots”  needed 
and  the  corresponding  computer  time  by  about  a  factor  of  five. 

The  second  step  involves  taking  the  above  data  and  constructing  the  reduced  order  model  which 
requires  about  1 5  minutes  of  computer  time. 

With  the  information  from  steps  one  and  two  above,  in  the  third  step  each  aeroelastic  solution  for 
a  given  Mach  number  and  set  of  structural  modes  and  one  combination  of  all  parameters  takes  a 
fraction  of  a  second.  This  provides  the  true  frequency  and  damping  (i.e.,  this  is  a  p  method)  for 
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each  aeroelastic  mode.  To  construct  a  root  locus  using  say  100  dynamic  pressure  values  takes 
less  than  one  minute. 

Flutter  Results 


Fig  4.2  shows  the  eigenvalue  root-loci  when  sweeping  through  various  mass  ratios  (to  which 
there  is  a  corresponding  flutter  speed  index)  when  solving  the  aeroelastic  eigenvalue  problem 
posed  by  the  reduced-order  aeroelastic  model  (Eq  (4.7))  for  various  Mach  numbers.  Solution 
snapshots  are  computed  for  the  first  five  given  wing  mode  shapes  for  reduced  frequencies 
(k  =  cob/ Ux)  from  k  =  0.0  to  k  =  0.5  in  Ak  =  0.1  increments.  This  configuration  flutters  for 
frequencies  less  than  0.5,  and  as  such,  solution  snapshots  for  k  >  0.5  are  unnecessary.  This 
results  in  a  total  of  55  available  POD  shape  vectors.  In  Fig  4.2,  the  curves  represent  the 
eigenvalues  corresponding  to  the  primarily  structural  natural  modes  as  mass  ratio  is  varied.  Our 
method  also  determines  the  aeroelastic  modes  originating  from  the  fluid  dynamic  modes  of  the 
POD/ROM.  For  the  range  of  mass  ratios  (0  <  p  <  500)  swept  through  in  these  parametric 
analyses,  the  fluid  dynamic  modes  are  very  damped,  and  as  such  lie  to  the  left  and  outside  of  the 
eigenspectrum  range  we  show.  As  can  been,  for  each  of  the  Mach  numbers,  the  first  structural 
mode  tends  to  be  the  critical  flutter  mode.  For  the  highest  Mach  number  however,  the  third 
structural  mode  can  go  unstable  if  the  mass  ratio  is  large  enough.  Also  from  this  figure,  it  is 
evident  that  it  is  unnecessary  to  use  all  55  of  the  available  POD  shapes.  If  fact,  with  less  than 
one  half  of  the  POD  modes  (25  for  instance),  relatively  converged  results  (in  the  sense  of  POD 
mode  refinement)  can  be  achieved. 

Fig  4.3  shows  the  computed  POD/ROM  flutter  speed  and  flutter  frequency  ratios,  along  with 
experimental  data  (Ref  38),  and  data  from  two  other  computational  methods  (Refs  40,  41),  as  a 
function  of  Mach  number. 

As  can  be  seen,  using  our  methodology,  we  produce  the  well  known  transonic  flutter  speed  dip, 
and  our  results  are  all  within  the  same  tolerance  to  the  experimental  results  as  the  other 
computational  methods.  Gupta  (Ref  41),  does  show  better  agreement  with  experiment  at  the  two 
supersonic  Mach  numbers,  and  Gupta  attributes  this  better  agreement  to  better  CFD  grid 
refinement.  In  future  work,  we  will  also  address  this  issue. 
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The  Use  of  Alternate  Modal  Excitations  for  Snapshots 


One  of  the  key  concerns  towards  in  the  POD/ROM  method  to  three-dimensional  flows  has  been 
whether  or  not  an  entire  set  of  solution  snapshots  must  be  computed  for  each  possible  structural 
configuration  of  interest.  That  is,  say  we  wish  to  consider  a  similarly  shaped  wing  that  has  a 
slightly  different  structural  configuration,  which  in  turn  means  the  wing  vibratory  mode  shapes 
are  different.  Does  this  mean  that  one  has  to  go  through  and  compute  a  whole  new  ensemble  of 
solution  snapshots  based  on  these  new  structural  motions  in  order  to  do  a  flutter  analysis  for  the 
new  wing  configuration.  Fortunately,  although  a  few  snapshots  based  on  the  new  modal  motion 
will  need  to  be  computed,  the  larger  number  of  snapshots  computed  at  numerous  frequencies  will 
be  unnecessary.  The  snapshots  computed  from  a  previous  wing  structural  configuration  will  still 
serve  the  purpose.  That  is,  a  couple  of  solution  snapshots  will  be  needed  for  the  new  structural 
motions,  however,  these  will  only  need  to  be  computed  at  the  end  points  of  the  frequency  range 
of  interest.  Fig  4.4  demonstrates  how  this  works. 

Fig  4.4  shows  the  real  and  imaginary  parts  of  the  coefficient  of  the  generalized  aerodynamic 
force  corresponding  to  the  first  mode  pressure  acting  through  the  first  mode  shape  as  a  function 
of  reduced  frequency  at  a  Mach  number  of  0.960.  The  coefficient  of  the  generalized 
aerodynamic  force  is  defined  as 


(4.8) 


where  qx  =  pmUl  /2  is  the  freestream  dynamic  pressure  and  cr  is  the  root  chord  length.  In  this 
definition,  pfo)  represents  the  frequency  dependent  unsteady  pressure  resulting  from  a  wing 
deformation  motion  of 


-  =  *,  (4.9) 

c. 

The  curves  presented  in  Fig  4.4  are  based  on  the  actual  solution  snapshots  and  thus  are  what  we 
desire  the  POD/ROM  to  model.  In  Fig  4.4.b,  the  POD/ROM  of  C0(  based  on  snapshots  for 
each  of  the  five  structural  mode  shapes  at  frequencies  of  £=0.0  and  £=1.0  (for  a  total  of  10 
snapshots)  is  compared  against  C&  _  for  the  actual  snapshots  for  all  frequencies  between  £=0.0 

and  £=1.0.  As  can  be  seen,  the  POD/ROM  matches  at  the  end  points  of  the  frequency  range  as  is 
expected,  however  this  crude  POD/ROM  performs  rather  poorly  for  the  intermediate 
frequencies.  Of  course,  if  we  use  snapshots  at  all  the  frequencies  between  £=0.0  and  £=1.0,  the 
POD/ROM  would  exactly  reproduce  the  data. 
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(a)  Snapshots 


(b)  POD/ROM  based  on  k=0.0  and  k=1.0 
modal  snapshots 


(c)  POD/ROM  based  on  basic  motions 
snapshots 


(d)  POD/ROM  based  on  k=0.0  and  k=1.0 
modal  snapshots  and  basic  motions 
snapshots 


Figure  4.4  Generalized  Force  Modeling  with  Unrelated  Mode  Shape  Snapshots 


Next,  in  Fig  4.4.C,  a  new  POD/ROM  for  Cfl  now  based  on  solution  snapshots  unrelated  to  the 

actual  mode  shapes  is  shown.  The  simple  wing  motion  snapshots  are  for  a  full  wing  plunge 
motion  (up/down),  full  wing  pitch  about  the  quarter  chord,  a  first  bending  type  of  motion  (wing 
is  fixed  at  the  root,  and  the  z  coordinate  component  of  deflection  varies  linearly  with  span),  and  a 
first  twist  type  of  motion  (wing  is  fixed  at  the  root,  and  the  pitch  varies  linearly  with  span)  for 
frequencies  from  &=0.0  to  A=1.0  at  Ak  =  0.1  increments  for  a  total  of  44  solution  snapshots.  As 
can  be  seen,  the  POD/ROM  in  this  case  also  perform  very  poorly.  Unbeknownst  however,  these 
solutions  are  in  fact  helping  to  reveal  the  dynamics  of  the  system.  In  fact,  when  one  uses  these 
snapshots  in  combination  with  the  actual  structural  mode  snapshots  solely  at  the  end  points  of  the 
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frequency  range  of  interest,  one  gets  a  POD/ROM  which  produces  very  accurate  results  to  ( 
as  is  evident  from  Fig  4.4.d. 

A  comparison  of  the  POD/ROM  Mach  number  flutter  trends  in  the  case  where  first,  the 
POD/ROM  is  based  on  solution  snapshots  corresponding  to  the  actual  modal  shapes  of  the  wing 
to  the  case  where  second,  the  POD/ROM  is  based  on  snapshots  using  the  simple  wing  motion 
mode  shapes  as  discussed  in  the  previous  paragraph  has  been  made.  The  results  are  virtually 
identical  at  all  but  the  highest  Mach  number.  There  is  some  difference  at  the  highest  Mach 
number,  again  suggesting  supersonic  flow  is  more  sensitive  for  this  wing. 

For  additional  details  of  the  work  presented  in  this  Chapter,  please  see  the  paper  (Ref  1 6) 
included  in  Appendix  C. 
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CHAPTER  5 


SIMULATION  OF  TRANSONIC  LIMIT  CYCLE  OSCILLATIONS  USING 

A  CFD  TIME-MARCHING  METHOD 


Summary 

We  have  conducted  a  transonic  LCO  investigation  for  an  NLR  7301  airfoil  using  NASA/Langley 
CFL3D/N-S  time-marching  approach  to  validate  with  measured  LCO/flutter  data  of  DLR.  The 
present  LCO  study  is  a  parallel  effort  to  that  of  Duke  University  whose  CFD  method  is  based  on 
a  frequency-domain  POD/ROM  EigenMode  approach. 

The  primary  objective  of  the  CFL3D/N-S  effort  was  aimed  at  to  investigate  the  impact  of 
turbulent  models  and  computing  time  on  a  viable  CFD  tool  for  transonic  LCO.  We  found  that 
the  LCO  solution  is  turbulent-model  dependent  and  time-step  size  sensitive.  These  issues  along 
with  the  possible  influence  on  LCO  due  to  initial  conditions  and  transition  points  remain  to  be 
subjects  for  further  research. 

In  using  the  Spalart-Allmaras  turbulent  model,  the  present  study  for  NLR7301  shows  that  the 
time-marching  CFL3D  computation  can  predict  LCO  and  with  good  frequency  agreement,  which 
to  a  large  extent  validates  the  transonic  test  result  of  Schewe  et  al  at  DLR.  However,  in  order  to 
obtain  a  fully  developed  LCO,  this  time-marching  computation  typically  requires  about  96  hours 
of  computer  time  on  a  1  GHz  computer.  In  view  of  its  long  computing  time,  a  CFD  time¬ 
marching  scheme  such  as  CFL3D/N-S  would  be  computationally  inefficient  as  a  tool  for  3D 
transonic  LCO  investigation.  The  above  observation  suggests  Duke’s  frequency-domain 
POD/ROM  EigenMode.  approach,  for  given  computer  limitation,  may  very  well  be  the  only 
viable  tool  for  3-D  transonic  LCO  investigation  at  present. 


Introduction 


Limit  Cycle  Oscillation  (LCO)  has  been  a  persistent  problem  on  several  current  fighter  aircraft 
and  is  generally  encountered  with  external  store  configurations.  Denegri  (Ref  1)  provided  a 
detailed  description  of  the  aircraft/store  LCO  phenomenon.  Norton  (Ref  2)  gave  an  excellent 
overview  of  LCO  of  fighter  aircraft  carrying  external  stores  and  its  sensitivity  to  store  carriage 
configuration  and  mass  properties. 

LCO  can  be  characterized  as  sustained  periodic  oscillations  which  neither  increase  or  decrease  in 
amplitude  over  time  for  a  given  flight  condition.  Using  an  s-domain  unsteady  aerodynamic 
model  of  the  aircraft  and  stores,  Chen,  Sarhaddi  and  Liu  (Ref  7)  has  shown  that  wing/store  LCO 
is  a  post-flutter  phenomenon  whenever  the  flutter  mode  contains  low  unstable  damping.  This 
type  of  flutter  mode  is  called  the  “ hump  mode".  Since  the  aircraft  structure  usually  contains 
structural  nonlinearity  such  as  friction  damping,  this  amplitude-dependent  friction  damping  can 
suppress  the  growth  of  amplitude,  thus  resulting  in  a  steady  state  oscillation.  This  is  known  as 
the  nonlinear  structural  damping  (NSD)  model  of  the  wing/store  LCO.  Although  not  thoroughly 
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proven  through  tests  or  numerical  simulations,  results  of  the  NSD  show  excellent  correlation 
with  flight  test  LCO  data  of  F-16  throughout  subsonic  and  transonic  Mach  numbers. 

On  the  other  hand,  other  long-standing  researchers,  notably  Meijer  and  Cunningham  (Ref  8), 
believe  that  the  wing/store  LCO  is  due  largely  to  the  transonic  shock  oscillation  and  shock 
induced  flow  separation,  called  Transonic  Shock/Separation  (TSS)  model.  Edwards  considers 
the  TSS  model  and  viscous  effects  are  two  major  factors  that  cause  transonic  LCO  for  wings.  He 
also  delineates  the  shock  buffet  phenomenon  as  opposed  to  that  of  transonic  LCO  (Ref  9).  It 
should  be  noted  that,  however,  there  is  no  conflict  in  the  NSD  model  and  the  TSS  model  in  that 
TSS  is  only  a  transonic  subset  of  NSD.  That  is  to  say  that  NSD  model  could  consistently 
interpret  the  LCOs  occurrence  at  subsonic  and  supersonic  speeds  as  well,  whereas  TSS  can  not. 

Recent  renewed  interest  in  LCO  is  perhaps  motivated  by  the  need  of  further  resolving  fighter 
LCO  and  the  current  advent  of  CFD  methodology  in  aeroelasticity.  There  are  two  potential 
computational  methods  for  LCO  prediction/investigation:  the  CFL3D.AE  code  (Ref  10) 
developed  and  supported  by  NASA/Langley  and  the  POD/ROM  EigenMode  approach  (Ref  11) 
originated  by  Dowell  and  Hall  of  Duke  University.  The  former  is  a  conventional  time-domain 
CFD  method  whereas  the  latter  a  frequency-domain  CFD  method,  using  aerodynamic  eigen 
modes. 

Objectives  o£CFD  Simulation 

Using  the  Navier-Stokes  option  of  CFL3D.AE  for  LCO  investigation,  ZONA  selects  an 
NLR7301  airfoil  (shock-free  design  condition  at  M  =  0.747,  C,  =  0.455)  as  the  case  studied. 
Primarily,  a  thorough  wind-tunnel  transonic  LCO/flutter  study  was  performed  for  a  2D 
supecritical  wing  with  NLR7301  airfoil  section  by  Schewe  and  his  associates  at  DLR  (Refs  42- 
44).  To  our  best  knowledge,  Schewe’s  work  is  perhaps  the  only  experimental  work  available  for 
measuring  2-D  transonic  LCO.  With  Schewe’s  data,  we  began  to  validate  our  CFD  simulated 
results  with  the  following  objectives: 

•  Attempt  to  understand  the  LCO  physics 

•  Investigate  the  impact  on  LCO  solutions  due  to  different  various  flow  parameters 

•  Evaluate  necessary  computer  resource  including  computing  time  required 

•  Establish  benchmarks  and  identify  an  efficient  3-D  CFD  method  for  LCO/flutter  prediction. 

LCO  Study  o£NLR  7301  Airfoil  Using  CFL3D 

Schewe  and  Deyhle  (Ref  42)  conducted  an  experiment  on  transonic  flutter  of  a  2-D  supercritical 
wing  with  an  NLR  7301  airfoil  section.  The  emphasis  of  this  experiment  was  to  investigate  the 
effect  of  the  flow  nonlinearity  on  the  aeroelastic  response  including  transonic  dip  and  transonic 
LCO.  From  the  aeroelastic  point  of  view,  an  aeroelastic  experiment  of  a  2-D  wing  is  more 
difficult  and  complicated  than  that  of  a  3-D  wing.  To  validate  a  CFD  methodology  with  such 
valuable  data  will  help  accomplishing  the  present  objectives.  The  following  is  an  account  of  our 
theoretical  modeling  of  the  experimental  set  up. 
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Test  Set-up  and  Equations  of  Motion 


Figure  5.1  Two-Degree-of-Freedom  Dynamic  Model 


Fig  5.1  depicts  a  simplified  model  of  the  2-degree-of-freedom  test  set-up.  The  2-D  wing  has  a 
chord  length  of  0.3  m  (c  =  0.3  m)  and  a  span  of  1  m  (b  =  1  m).  The  pitching  spring  and  heaving 
spring  are  attached  to  the  same  c/4  position.  The  corresponding  2-degree-of-freedom  equation  of 
motion  of  the  set-up  reads, 

mh  ~sa  1  f  h  \  Dh  0  I  [  1  Kh  0  I  f  1  f  L(t)  1 

i  r +  i  r +  i  r =  1  r  (5-D 

_-5a  /c/4  J  l  «  J  L°  Da  ][d\  L0  Ka\{ct)  M{t)  J 

where 

mh  is  the  total  mass  ( mh  =  26.64  kg) 

Ir,A  is  the  mass  moment  of  inertia  about  c/4  ( Ic ,  4  =  0.086  kg.m2) 

su  is  the  static  unbalance  (sa  =  0.378  kg.m) 

Dh  and  Da  are  the  damping  factors  of  the  heave  motion  (h)  and  the  pitch  motion 
(a),  respectively  (Dh  =  82.9  kg/s  and  Da  =0.197  kg.m2/(rad-s)) 

Kh  and  Ka  are  the  stiffness  of  the  heaving  spring  and  pitching  spring,  respectively 
(Kh  -  1.21  x  106  N/m  and  Ka  =  6.68  x  103  N-m/rad),  and 

L(t)  and  M(t)  are  the  aerodynamic  lift  and  moment,  respectively. 

The  numerical  values  of  the  structural  terms  in  Eq  (5.1)  can  also  be  found  in  Ref  43.  To  perform 
the  time-marching  CFD  computation,  it  is  convenient  to  convert  Eq  (5.1)  from  the  physical 
degrees  of  freedom  to  the  modal  coordinates,  i.e.: 

•*[  =  [«>]{«}  (2) 
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where  q  is  modal  coordinate  and  <j>  is  the  modal  matrix  of  the  undamped  structure 


1 *  = 


-0.1735  0.1004 

0.9277  3.403 


Substituting  Eq  (5.2)  into  Eq  (5.1)  and  pre-multiplying  the  resulting  equation  by  <pJ  yield: 


W  tel  + 


'2  (ohCh 

0 

fc}+ 

6)h  2  0  " 

W  =  q*bf< 

cCL(t) 

0 

2  <OaCa_ 

i 

o 

K 

L _ 

1 

(5.3) 


where 

coh  and  coa  are  the  undamped  natural  frequencies  of  the  heaving  and  pitching 
motions,  respectively  ( coh  =205.4  rad/s  and  coa  =299.5  rad/s) 


Ch  and  C,a  are  the  heaving  and  pitching  damping  ratios,  respectively  (t,h=  0.00648 

and  Qa  =  0.00474).  Note  that  the  off-diagonal  terms  in  the  damping 
matrix  are  assumed  to  be  zero  for  simplicity. 

qoo  is  the  dynamic  pressure,  and 

Ci(t)  and  Cm(t)  are  the  non-dimensional  aerodynamic  lift  and  moment  coefficients  of  the 
2-D  section  of  the  wing,  i.e.: 


CL(t)  = 


m 

be 


and 


Cm{t)  = 


M(f) 

q*bc2 


(5.4) 


Ci(t)  and  Cm(t)  are  computed  by  the  CFL3D  code  based  on  the  NLR 
7301  full  chord  length. 


•  Steady  State  Results 

Fig  5.2  shows  a  C-type  grid  with  273x93  mesh  points  around  the  NLR  7301  airfoil.  For 
validation,  we  first  performed  a  2-D  steady  computation  with  the  Baldwin-Lomax  and  Spalart- 
Allmaras  turbulence  models  assuming  fully  turbulent  flow 
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Figure  5.2  C-Type  Grid  Around  NLR  7301  Airfoil  (273  x  93) 

The  comparison  of  the  computed  pressure  coefficient  (Cp)  with  the  experimental  data  is  shown  in 
Fig  5.3.  Good  agreement  between  these  two  sets  of  CP  distribution  can  be  seen  except  at  the 
20%  chord  on  the  suction  side.  Weber  et  al  (Ref  45)  presented  a  similar  result  and  suggested  that 
this  fully  turbulent  result  can  be  improved  by  taking  a  fixed  transition  into  account.  However, 
due  to  the  lack  of  test  data  in  the  transition  onset  location,  no  fixed  transition  model  is  considered 
in  the  present  steady  and  unsteady  aerodynamic  computations. 


Figure  5.3  Steady  Pressure  Distribution 
CFL3D:  Mach  No.=  0.753,  AOA=-0.08°,  Re=l. 727x1 06 
Experiment:  Mach  No.=  0.768,  AOA=1.28°,  Re=l. 727x1 06 
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Time  Step  Size  for  Convergence 


It  is  well  known  that  the  transonic  flow  characteristics  of  the  supercritical  airfoil  is  very  sensitive 
to  off-design  conditions.  The  position  and  strength  of  the  transonic  shock  can  change  rapidly 
due  to  a  small  angle-of-attack  and/or  Mach  number  perturbation  (Ref  46).  This  indicates  that  for 
transonic  unsteady  flow  computation  of  supercritical  airfoils,  the  time  step  must  be  kept 
sufficiently  small  and  the  number  of  Newton  sub-iterations  within  each  time  step  must  be 
adequate  to  ensure  the  solution  convergence.  On  the  other  hand,  this  stringent  condition  for  low 
speed  unsteady  flow  computation  can  be  largely  relaxed. 

To  show  this,  we  conducted  a  time-step  size  study  for  solution  convergence  on  the  NLR  7301 
airfoil  at  M=0.05  and  M=0.753  by  imposing  a  sinusoidal  motion  with  oscillatory  frequency  at  50 
Hz  and  pitching  amplitude  of  1  degree. 

For  the  transonic  case  (M=0.753),  the  time  histories  of  the  lift  coefficient  with  N=200/subit=3, 
N=800/subit=6,  and  N=800/subit=9  is  shown  in  Fig  5.4.a,  where  N  is  the  number  of  time  steps 
within  each  cycle  and  “subit”  is  the  number  of  Newton  subiterations  within  each  time  step.  It 
can  be  seen  that  the  solution  varies  largely  from  N=200/subit=3  to  N=800/subit=6  and  thereupon 
it  varies  little  up  to  N=800/subit=9,  suggesting  that  a  time  step  size  corresponding  to 
N=800/subit=9  must  be  adopted  for  transonic  unsteady  aerodynamic  computation. 

The  same  type  of  time  histories  for  the  low  speed  unsteady  flow  (M=0.05)  is  shown  in  Fig  5.4.b, 
but  with  N=800/subit=9,  N=200/subit=3,  and  N=25/subit=3.  Little  variations  of  the  solutions 
from  N=800/subit=9  to  N=25/subit=3  can  be  seen,  suggesting  that  a  time-step  corresponding  to 
N=25/subit=3  is  sufficient  for  low  speed  unsteady  flow  computation,  this  also  indicates  that  the 
transonic  LCO  analysis  requires  at  least  two  orders  of  magnitude  more  computing  time  than  that 
of  a  low  speed  analysis.  In  fact,  in  the  following  sections  we  will  show  that  to  obtain  a  fully 
developed  transonic  LCO  time  history  takes  several  days  of  computer  time  on  a  1  GHZ 
computer. 
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Number  of  Iterations 


(a)  M=0.753 


Nondimensional  Time 
(b)  M=0.05 


Figure  5.4  Time  Step  Size  for  Solution  Convergence 


•  Steady  AoA  Condition  and  Static  Aeroelastic  Equilibrium  (SAE)  Condition 

Based  on  the  complete  Eq  (5.3),  the  time-marching  computation  can  then  follow  with  the 
resulting  SAE  condition  being  the  initial  condition.  Here,  the  initial  conditions  in  terms  of 
generalized  coordinates  are  imposed  at  the  values  of  qi=  0  and  q2  =  q2(0)/2  ,  see  Table  5.2.  In 
passing,  we  note  that  different  initial  conditions  could  affect  the  resulting  LCO  solution  in  the 
presence  of  flow  nonlinearity.  Nonetheless,  we  merely  leave  the  question  of  the  initial- 
condition  influence  on  LCO  as  a  future  research  topic. 
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In  Schewe’s  experiment,  there  exist  two  starting  flow  conditions:  namely,  a  steady  AoA 
condition  and  a  static  aeroelastic  equilibrium  (SAE)  condition.  A  steady  AOA  condition  is 
defined  by  the  initially  selected  flow  condition,  whereas  the  SAE  condition  is  the  subsequent 
measured  mean  AoA  position  during  the  steady  state  oscillation.  In  terms  of  our  numerical 
simulation,  the  SAE  condition  is  one  which  obtained  from  solving  the  static  aeroelastic  version 
of  Eq  (5.3),  i.e.,  by  imposing  a  very  large  structural  damping  in  order  to  nullify  the  velocity  and 
acceleration  terms. 

Presented  in  Table  5.1  are  the  starting  flow  conditions  of  Schewe  and  those  of  two  simulated 
methods.  The  difference  in  the  test  and  simulated  values  of  the  steady  AoA  condition  results 
from  correlation  of  the  corresponding  Cp’s  and  lift  coefficients.  The  slight  discrepancy  in  the 
SAE  AoA  between  two  simulated  approaches  remained  to  be  clarified.  Table  5.2  presents  the 
SAE  generalized  coordinates  as  converted  from  the  resulting  SAE  flow  conditions. 


Table  5.1  Two  Starting  Flow  Conditions 


Steady 

AOA 

Condition 

Static  Aeroelastic 

Equilibrium 
(SAE)  Condition 

Knipfer  et  al  (Ref  44) 

M  =  0.768,  a  =  1.91° 

M  =  0.768,  a  =  1.28° 

Weber  et  al  (Ref  45) 

M  =  0.753,  a  =  0.6° 

M  =  0.753,  a  =  -  0.08° 

Present 

M  =  0.753,  a  =  0.6° 

M  =  0.753,  a  =  0.078° 

Table  5.2  SAE  Generalized  Coordinates 


Dynamic  Pressure 

Generalized  Coordinates 

qi(0) 

q2(0) 

Pdynamic  “15  kPa 

-  0.0059 

-0.0031 

^dynamic  “  12.6  kPa 

-0.0051 

-  0.0023 

Pdynamic  “  9.5  kPa 

-  0.0045 

-0.0017 
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•  Verification  of  Flutter  Boundary  at  M  —  0. 753 

Fig  5.5  presents  Schewe  and  Deyhle’s  experimental  flutter  boundary  for  a  2-D  supercritical  wing 
with  NLR  7301  airfoil  section  (Ref  42).  In  the  subsonic  range  up  to  M  =  0.74,  the  critical  value 
of  dynamic  pressure  is  nearly  constant.  Transonic  dip  occurs  when  M  >  0.74  where  a  dramatic 
drop  of  the  critical  dynamic  pressure  can  be  seen.  In  this  transonic  dip  region,  Schewe  and 
Deyhle  reported  that  various  types  of  LCOs  were  observed.  Using  a  Fourier  transformation  of 
the  Navier-Stokes  solver  generated  unsteady  aerodynamics,  Knipfer  and  Schewe  (Ref  43) 
performed  a  frequency-domain  flutter  analysis  and  concluded  that  these  LCO  were  genuine 
flutter  cases.  This  conclusion  can  be  further  supported  by  the  fact  that  the  LCO  appears  already 
at  Mach  number  below  the  buffet  limit.  To  verify  this,  we  performed  two  CFL3D  aeroelastic 
computations  at  M  =  0.753,  one  at  q,*,  =  9.5  kPa  and  the  second  one  at  q*  =  15  kPa.  In  both 
computations,  the  speed  of  sound  is  fixed  at  254.7  m/s  and  invariant  to  the  air  density. 

The  time  histories  of  the  heaving  motion  (h)  and  the  pitching  motion  (a)  of  these  two  cases  are 
presented  in  Figs  5.6  and  5.7,  respectively.  As  expected,  divergent  motions  at  q*,  =  15  kPa  and 
convergent  motions  at  q»  =  9.5  kPa  are  obtained.  This  confirms  that  the  LCO  of  the  2-D 
supercritical  wing  with  NLR  7301  airfoil  section  is  indeed  a  post-flutter  phenomenon,  which 
initiates  from  a  classical  flutter  instability. 


Mach  Number 


Figure  5.5  Flutter  Boundary:  Dynamic  Pressure  vs.  Mach  Number 
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Figure  5.6  Time  History  at  Mach  0.753,  AOA=0.6°,  q=9.5  kPa,  Re=1.727xl06 
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0  0.5  1  1.5 

Time  (seconds) 

Figure  5.7  Time  History  at  Mach  0.753,  AOA=0.6°,  q=15  kPa,  Re=1.727xl06 


•  LCO  Study  at  M  =  0. 753  with  Various  Turbulence  Models 

LCO  is  generally  caused  by  the  nonlinearity  of  a  dynamic  system.  The  structural  set  up  of  the  2- 
D  supercritical  wing  section  is  basically  linear,  whereas  the  LCO  observed  in  the  experiment 
results  from  unsteady  nonlinear  transonic  flow  including  effects  of  transonic  shock  movement 
and  possibly  shock-induced  separation.  This  suggests  that,  in  order  to  correlate  with  the 
experimental  LCO  results,  the  impact  of  different  turbulence  models  on  the  predicted  LCO  must 
be  studied  first.  To  do  this,  we  first  performed  a  CFL3D  computation  with  the  Baldwin-Lomax 
model  at  M  =  0.753,  q*.  =  12.6  kPa  and  Re  =  1.727  x  106.  Since  the  qoo  =  12.6  kPa  condition  is 
already  beyond  the  flutter  boundary,  an  initially  divergent  motion  at  small  amplitude  is  expected 
from  the  CFL3D  result.  This  is  verified  by  the  resulting  time  histories  of  the  h  and  a  motions 
presented  in  Fig  5.8  where  the  divergent  motions  up  to  t  =  1.0  sec  can  be  clearly  seen.  At  t  =  1.2 
sec,  a  sudden  increase  of  amplitude  occurs  but  the  divergent  rate  of  the  motion  decreases. 
However,  when  further  extending  the  time-marching  computation  from  t  =  1.5  to  2.8  sec,  the 
divergent  motion  slows  down  yet  still  sustaining  and  does  not  reach  a  steady  state  oscillation, 
after  96  hours  computation  on  a  1  GHz  computer.  In  many  instances  of  such  LCO  computing 
tasks,  the  difficulty  we  faced  is  that  a  once  started  the  computation  must  be  carried  on  typically 
for  several  days  until  a  LCO  reveals  itself.  Otherwise,  the  computation  shall  be  terminated  by 
computing  time  limitation,  if  a  LCO  has  not  yet  been  reached  by  then;  this  undetermined  LCO  is 
marked  as  a  divergent  case. 
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Next,  the  CFL3D  computation  was  repeated  at  the  same  flow  condition  but  using  the  Spalart- 
Allmaras  turbulence  model  and  the  result  is  shown  in  Fig  5.9.  Similar  type  of  divergent  motion 
up  to  t  =  1.0  sec  to  that  of  the  Baldwin-Lomax  model  is  obtained.  But  this  time,  a  steady  state 
oscillation  is  finally  reached  by  extending  the  time-marching  computation  up  to  t  =  2.4  sec.  The 
predicted  LCO  frequency  is  approximately  34  Hz  which  correlates  well  with  that  of  the 
experiment.  But  the  predicted  LCO  pitch  amplitude  is  approximately  4°  which  is  off  by  an  order 
of  magnitude  comparing  to  the  experimental  LCO  amplitude  (what  Schewe  measured  is  0.6°). 
Similar  disagreement  in  terms  of  the  LCO  amplitude  was  found  by  Weber  et  al  (Ref  45).  Closer 
agreement  of  the  LCO  amplitude  was  obtained  by  Castro  et  al  (Ref  47)  where  the  wind  tunnel 
interference  effects  was  included  in  their  CFD  simulation.  Note  that  the  interference  effect  is  not 
considered  in  the  present  LCO  study  and  that  of  Weber  et  al. 
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Figure  5.8  LCO  Study  Using  Baldwin-Lomax  Model 
at  Mach  0.753,  q=12.6  kPa,  Re=1.727xl06. 
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CHAPTER  6 


PHASES  II  /  III  PLAN 

l 

Summary 

The  ZONA  Technology/Duke  University  team  proposes  a  challenging  set  of  tasks  to  be  accomplished  in  two  years 
time  that  builds  upon  the  very  substantial  progress  already  made  in  Phase  I.  Specifically,  the  2D,  Euler  Equations 
CFD  model  (developed  in  Phase  I)  that  uses  a  harmonic  balance  solution  technique  will  be  extended  to  include  the 
effects  of  viscosity  (Navier-Stokes  equations)  and  multiple  structural  modes.  In  parallel  the  3D,  Euler  Equations 
CFD  time  linearized  model  (also  developed  in  Phase  I)  which  already  includes  multiple  structural  modes  and  uses  a 
POD/ROM  solution  method  will  be  extended  to  now  include  full  dynamic  nonlinear  effects  using  a  harmonic 
balance  solution;  and  then  further  extended  to  include  the  effects  of  multi-bodies  (e.g.  wings  plus  tip  missiles  and 
stores).  With  each  new  advance  in  modeling  and  solution  technique,  flutter  and  LCO  analysis  and  prediction  will  be 
made  and  the  results  compared  to  existing  data  for  2D  flow  over  airfoils  and  3D  flow  over  wings  and  wing/stores. 
The  culminating  prediction  and  analysis  will  be  for  the  F-16  and  F/A-18  aircraft.  Based  upon  the  very  encouraging 
results  of  Phase  I  using  the  POD/ROM  and  Harmonic  Balance  solution  techniques  and  in  comparison  with 
benchmarking  calculations  using  a  state-of-the-art  time  marching  3D  flow  solver  (CFL3D),  we  anticipate  that  the 
computer  models  to  be  developed  in  Phase  II  (like  their  Phase  I  predecessors)  will  enjoy  a  very  substantial 
computational  efficiency  advantage  over  other  existing  unsteady  CFD  models.  Typical  computer  time  reductions  are 
expected  to  be  a  least  two  orders  of  magnitude  for  LCO  calculations  and  even  more  for  flutter  boundary  analysis 
per  se. 


ZONA  Technology,  Inc.  and  Duke  University,  the  ZONA  team,  jointly  hereby  propose  a  three- 
phase  global  program  entitled  “ Nonlinear  Reduced-Order  Modeling  of  Limit  Cycle  Oscillations 
of  Aircraft  Wings  and  Wing/Stores".  The  overall  objective  of  the  proposed  program  is  to  develop 
the  frequency-domain  POD/ROM/HB  EigenMode  method  in  order  to  further  the  understanding, 
accurate  and  efficient  prediction,  and  control  of  LCO/flutter  for  aircraft  wings  and  wing/stores. 

6.1  Program  Goals 

The  overarching  goal  for  the  Phase  II  work  is  to  build  on  the  very  substantial  progress  made  in 
Phase  I  and  to  develop  a  capability  for  a  highly  computationally  efficient  and  physically  accurate 
mathematical  modeling  of  limit  cycle  oscillations  (LCO)  and  other  nonlinear  aeroelastic 
phenomena.  The  approach  uses  the  concepts  of  aerodynamic  modes  as  well  as  structural  modes. 
Such  models  provide  substantially  improved  physical  understanding  and  also  more  accurate 
prediction  of  LCO  in  particular  and  nonlinear  aeroelastic  responses  in  general.  This  capability 
may  also  lead  to  a  rational  prediction  of  buffet  onset  due  to  aerodynamically  nonlinear 
oscillations.  Such  models  will  include  the  following  physical  effects: 

•  Three  dimensional  (3D)  as  well  as  2D  flows 

•  Nonlinear  as  well  as  time  linearized  aerodynamic  pressures  and  forces 

•  Multi-body,  e.g.,  wings  plus  stores,  as  well  as  single  body  (e.g.,  wing  only)  configurations 

•  Many  as  well  as  few  structural  modes 

•  LCO  (nonlinear)  as  well  as  flutter  (linear)  aeroelastic  predictions 

6.1.1  The  Roadmap:  Overall  Program  Approach 

In  Fig  6.1,  a  roadmap  to  achieve  this  ambitious  set  of  goals  is  shown.  Concluding  with  Phase  III, 
all  of  the  above  capability  will  be  achieved. 
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Figure  6.1  The  Roadmap:  Overall  Program  Approach 
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6.1.2  Brief  Description  of  the  Roadmap 

Phase  I  included: 

•  Feasibility  study  of  a  time-marching  CFD/N-S  method 

•  Development  of  the  frequency-domain  POD/ROM  EigenMode  and  Harmonic  Balance 
methods  including  2-D  nonlinear  LCO/flutter/ffeeplay  case  studies 

•  Development  of  the  3-D  time-linearized  frequency-domain  POD/ROM  EigenMode  method 
including  a  flutter  case  study 

Phase  II  will: 

•  Develop  the  2-D  N-S  version  of  the  above  proposed  approach 

•  Extend  the  geometry  capability  to  3-D  including  wing/store  configurations 

•  Generalize  the  2-D/3-D  methods  to  include  Harmonic  Balance  technique  for  nonlinear  LCO 
analysis 

•  Include  multiple  structural  mode  analysis  capability 

•  Perform  case  studies  including  F-16  and  F-l  8  wing/store  LCO  and  business  jet  wing  LCO 
Phase  III  will  be  funded  by  ZONA  with  the  following  goals: 

•  Accomplish  the  final  step  of  the  development  -  fully  develop  a  nonlinear  3-D  Navier-Stokes 
version  of  the  proposed  method  for  complex  geometry 

•  Perform  LCO  for  fighters  F-16  and  F-l 8  with  stores  including  interaction  with  structural 
nonlinearity 

•  Combine  computational  advantages  of  POD/ROM  and  Harmonic  Balance  methods 

•  Once  fully  developed,  a  suite  of  computer  codes  will  be  packaged  as  a  software  product, 
ready  for  commercialization 

6.1.3  Discussion  of  Phases  I,  II  and  III 

•  Phase  I 

In  Phase  I  we  have  put  in  place  two  essential  parts  of  this  capability.  One  is  that  we  have 
developed  a  time-linearized  3D  aerodynamic  model  based  upon  a  modal  representation  using  the 
method  of  Proper  Orthogonal  Decomposition(POD).  With  this  aerodynamic  model  we  have 
performed  a  transonic  flutter  analysis  for  the  AGARD  445.6  wing  and  compared  these  results  to 
those  available  in  the  literature.  Our  results  show  good  correlation  with  data  obtained  by 
previous  investigators  while  reducing  the  computational  cost  for  the  aeroelastic  analysis  by 
several  orders  of  magnitude  over  more  conventional  time  marching  CFD  analyses.  In  this 
method,  an  accurate  representation  of  the  true  aeroelastic  eigenvalues  and  eigenvectors  is 
obtained  and  thus  physically  accurate  root  loci  may  be  obtained.  This  aeroelastic  eigenmode 
information  is  very  valuable  for  both  improved  physical  understanding  and  for  use  by  engineers 
who  will  be  designing  active  control  devices  for  aeroelastic  systems  with  smart  materials  or  more 
conventional  control  technology.  Typically  we  find  fewer  than  50  aerodynamic  states  are  needed 
and  often  many  fewer  states  may  be  used  with  acceptable  accuracy.  The  original  CFD  models 
from  which  the  aerodynamic  modes  are  determined  have  tens  of  thousands  of  states  in  2D  flow 
and  hundreds  of  thousands  or  more  in  3D  flow. 
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Also  in  Phase  I  we  have  developed  a  dynamically  nonlinear  aerodynamic  model  for  2D  flow 
based  upon  the  harmonic  balance  method  (a  Fourier  series  representing  the  multi-harmonics  of  a 
nonlinear  system).  This  aerodynamic  model  has  been  used  to  analyze  LCO  and  takes  into 
account  both  fluid  and  structural  nonlinearities.  These  studies  include  the  first  systematic  study 
of  LCO  due  to  control  surface  freeplay  in  the  transonic  regime  and  also  have  provided  new 
physical  insight  into  LCO  arising  from  aerodynamic  nonlinearities  due  to  large  shock  motion  for 
single  degree  of  freedom  flutter  as  well  as  aeroelastic  divergence.  Three  abstracts  have  been 
prepared  and  submitted  to  the  2001  AIAA  SDM  conference  describing  the  Phase  I  work  and 
these  are  included  here  as  appendices  for  the  convenience  of  the  reader. 

•  Phase  II 

The  challenge  now  is  to  construct  and  follow  a  roadmap  that  will  lead  to  the  ultimate  capability 
that  is  desired.  To  that  end,  consider  Fig  6.1  again  and  now  focus  on  the  proposed  Phase  II 
effort.  In  Year  1  of  Phase  II,  the  2D  nonlinear  aerodynamic  inviscid  flow  model  developed  in 
Phase  I  will  be  extended  to  include  viscous  effects  and  LCO  analyses  will  be  done  for  a  few 
structural  degrees  of  freedom.  This  will  lead  to  a  significant  advance  toward  our  goal  of 
including  viscous  flow,  nonlinear,  and  multi-structural  mode  capabilities.  In  the  same  year  our 
3D,  time  linearized  aerodynamic  model  will  be  extended  to  include  inviscid  nonlinear  effects. 
This  latter  work  also  builds  upon  the  achievements  of  Phase  I  where  3D  and  nonlinear  effects 
were  first  considered  separately.  In  general,  our  approach  to  the  road  map  is  to  add  a  new 
capability  (such  as  viscous  flow  effects  or  nonlinearities)  first  in  a  2D  model  and  then  to  use  this 
new  understanding  as  a  basis  for  adding  the  same  capability  to  the  3D  flow  models. 

In  year  2  of  Phase  II,  the  2D  aerodynamic  model  will  be  brought  together  with  a  multi-structural 
model  to  demonstrate  the  feasibility  of  doing  LCO  calculations  with  a  many  mode  model  for 
both  the  fluid  and  the  structure.  Harmonic  balance  analyses  are  often  used  with  a  relatively 
small  number  of  modes,  but  this  work  will  break  new  ground  for  models  of  the  size  typical  of 
aeroelastic  systems,  i.e.  10-50  structural  modes  and  a  comparable  number  of  aerodynamic 
modes.  In  the  same  year  the  3D  flow  model  will  be  extended  to  include  multi-bodies,  e.g.  wing 
plus  stores.  Such  an  extension  is  most  important  in  the  context  of  3D  flows. 

However  this  represents  a  larger  technical  risk  than  if  we  were  to  first  do  this  work  in  2D  flow. 
Nevertheless  we  are  confident  that  this  work  can  and  be  done.  LCO  analyses  will  be  performed 
with  this  code  once  it  is  available. 

•  Phase  III 

With  the  accomplishments  of  Phase  II,  years  1  and  2,  in  hand  a  strong  foundation  for  the  work  of 
Phase  III  will  be  laid.  Also  the  capability  already  developed  in  Phase  I  and  the  additional 
capability  to  be  developed  in  Phase  II  is  itself  quite  substantial  and  will  represent  a  marked 
advance  over  current  methods  for  LCO  and  other  aeroelastic  analyses.  In  the  following  sections 
are  discussed  two  major  themes  of  the  proposed  new  work.  The  first  theme  is  concerned  with 
novel  methods  for  determining  LCO  for  multi-mode  structural  systems  with  the  aerodynamic 
forces  represented  by  2D  or  3D  POD/ROM.  The  second  theme  addresses  the  question  of  how 
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we  plan  to  add  viscosity  to  our  POD/ROM.  With  respect  to  the  POD/ROM  representation  of 
more  complex  wing  and  body  geometries,  our  initial  approach  will  be  to  acquire  and  use  existing 
methodologies  for  conventional  CFD  codes,  e.g.,  the  needed  more  complex  computational  grids, 
and  then  simply  apply  our  POD/ROM  methodology  to  those  more  elaborately  gridded  CFD 
models.  We  anticipate  this  approach  will  be  successful.  However,  we  are  prepared  to  consider 
enhancements  of  existing  CFD  methods,  if  that  proves  to  be  needed.  For  the  determination  of 
LCO  and  the  inclusion  of  viscosity  we  anticipate  the  need  for  new  methodology  developments 
and  hence  those  are  discussed  in  section  5.0. 

6.1.4  The  ZONA/Duke  Team 

ZONA  Technology,  Inc.  (ZONA)  and  Duke  University  (Duke)  have  formed  a  strong  team  with  a 
comprehensive  background  to  handle  this  challenging  project.  Professor  Hall  (P.I.)  and 
Professor  Dowell,  one  of  the  world’s  leading  aeroelasticians,  at  Duke  are  the  originators  of  the 
proposed  ROM/POD  method,  which  has  gained  much  attention  in  the  aerospace  community  in 
recent  years. 

P.C.  Chen  (P.I.)  and  Danny  D.  Liu  at  ZONA  are  the  small  business  counterpart  for  this  project. 
For  over  a  decade,  ZONA  continues  to  perfect  the  proven  technology  in  aeroelastic  software 
such  as  the  ZONA51  code  in  MSC/NASTRAN  (Aero  Option  II),  now  an  industrial  standard  for 
supersonic  aeroelastic  analysis.  Meanwhile,  ZONA  has  an  outstanding  record  in  handling  U.S. 
government  contracts  as  evidenced  by  ZONA’s  successful  performance  in  five  ongoing  contracts 
(since  1995):  supported  by  AFRL/VA  (on  ASTROS*  and  VSS),  by  NAWC/Navy  (on  Missile 
Fin  Aeroelastic  Tailoring),  and  by  NASA/Langley  (on  BEM  Solver  for  CFD/CSD  Interfacing), 
by  NAVAIR/Navy  (on  Reconfigurable  Adaptive  Control  of  LCO  and  ASE  instability),  by 
NAVSEA/Navy  (on  ERGM  Projectile  Smart  Structure  Control).  ZONA  was  awarded  three 
SBIR/Phase  II  contracts  from  the  above  organizations  in  1999.  This  STTR  effort  between 
ZONA  and  Duke  is  supported  by  Boeing  as  evidenced  by  the  following  endorsement  letter  from 
Mr.  Rudy  Yurkovich. 

6.2  Phase  II  Objective  and  Anticipated  Benefits 

The  overall  Phase  II  technical  objective  is  to  develop  a  highly  computationally  efficient  and 
physically  accurate  mathematical  modeling  of  nonlinear  aeroelastic  phenomena.  Such  models 
provide  substantially  improved  physical  understanding  and  also  more  accurate  prediction  of 
nonlinear  aeroelastic  phenomena  such  as  LCO. 

Technical  Objectives^ 

Based  on  the  very  substantial  progress  made  in  Phase  I,  we  will  achieve  the  following  specific 
objectives  in  Phase  II: 

1 .  Development  of  a  2-D  viscous  version  of  the  harmonic  balance  technique. 

2.  Extension  of  the  geometry  capability  to  3-D  including  aircraft/store  configurations. 

3.  Generalization  of  the  2-D/3-D  methods  to  include  harmonic  balance  technique  for  LCO 
analysis  involving  aerodynamic  and  structural  nonlinearities. 
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4.  Inclusion  of  multiple  structural  mode  capability  for  large  finite  element  models  of  aircraft 
with  stores. 

5.  Validation  of  the  above  methodologies  by  performing  case  studies  including  F-16  and 
F/A- 18  wing/store  LCO  and  a  business  jet  wing  LCO. 

Anticipate^  Benefits 

Once  the  Phase  II  technical  objective  is  achieved,  its  anticipated  benefits  are: 

1 .  The  technical  advantages  of  the  proposed  approach  over  the  conventional  time-marching 
CFD  method  are  as  follows: 

•  Computational  efficiency 

•  Nonlinear  LCO  prediction  capability 

•  True  aeroelastic  damping  solutions  can  be  obtained 

2.  Once  fully  developed,  the  proposed  method  will  benefit  the  AF  and  the  aerospace 
community  in  the  following  areas; 

•  Becomes  an  efficient  aeroelastic/flutter  tool  for  the  industry 

•  Promotes  an  understanding  in  the  origins  of  LCO,  thus  providing  effective  means  to 
control  LCO 

•  Prediction  technique  will  help  resolve  the  wing/store  LCO  of  F-16  and  F-18  aircraft,  a 
solution  urgently  needed  by  Lockheed  Martin,  Boeing  and  AF. 

•  Address  the  transonic  tail  buzz  problem  of  many  next  generation  aerospace  vehicles  (e.g., 
F-22,  JSF,  X-33,  X-34,  etc.) 

•  Provides  a  sound  foundation  and  efficient  database  in  frequency-domain  for  effective 
transonic  aeroservoelasticity  (ASE)  and  MDO  applications. 

6.3  Phase  II  Work  Plan 

The  ZONA/Duke  team  proposes  a  two-year  effort  to  accomplish  the  technical  objectives  of 
Phase  II.  In  order  to  clearly  illustrate  the  Phase  II  work  plan  and  its  interrelationship  with  Phase 
I  and  Phase  III,  we  present  an  overall  program  task  chart  shown  in  Fig  6.2.  This  overall  program 
task  chart  is  derived  from  the  roadmap  of  Section  2.0  and  explicitly  specifies  the  tasks  involved 
in  the  overall  program 
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TL  =  Time-Linearized  (Euler),  HB  =  Harmonic  Balance  (Nonlinear), 

N-S  =  Navier-Stokes  (Add  viscosity),  SM  =  Several  Structural  Modes, 
MM  =  Multiple  Structural  Modes,  SW  =  Single  Wing,  W/S  =  Wing/Store 


Figure  6.2  Alternative  Chart  to  Roadmap 

The  Phase  II  development  is  planned  according  to  guidelines:  (see  Fig  6.1  Roadmap) 

•  Nonlinearities  that  produce  LCO  are  of  two  types:  aerodynamic  and  structural.  In  Phase 
II,  we  will  develop  models  that  can  analyze  both  types  of  nonlinearities  in  a  unified  way. 
At  the  end  of  Phase  II,  we  will  have  the  capability  to  analyze  cases  in  which  the  structure 
alone  is  nonlinear,  and  both  the  fluid  and  structure  are  nonlinear. 

•  The  complexity  of  the  aeroelastic  system  increases  as  Phase  II  advances: 

Geometry.  2-D  -»  3-D/Single  Wing  3-D  Wing/Store 
Equation  Hierarchy.  Euler  — >  Navier-Stokes 

-  Nonlinearity.  Time  Linearized  Approach  ->  Harmonic  Balance  Approach 

-  Mode :  Single  Degree  of  Freedom  -*  Multiple  Degree  of  Freedom 

-  Structural  Nonlinearity.  Linear  Structure  ->  Nonlinear  Structure 

In  this  section,  the  theoretical  formulation  of  the  methodologies  shown  in  Fig  6.1  will  be 
presented  in  details  (Sections  6.3.1  -  6.3.7).  Selected  test  cases  involving  a  single  wing  and  two 
fighter  aircraft  with  stores  configurations  are  presented  in  Section  6.3.8.  Phase  II  statement  of 
work  consisting  of  1 1  main  tasks  are  discussed  in  Section  6.3.9.  Program  Schedule  and 
Deliverables  are  shown  in  Sections  6.3.10  and  6.3.1 1,  respectively. 
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6.3.1  Three-Dimensional  Navier-Stokes  and  Euler  Equations 

We  start  with  the  top-down  approach: 

The  3-D  Reynolds-Averaged  Navier-Stokes  Equations  read: 

au  .  a(F-Fv)  ,  3(G-G„)  .  a(H-Hv)  _  „ 
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and  p,  p,  (u,  v,  w),  e,  h  are  the  flow  density,  pressure,  velocity,  energy  and  enthalpy, 
respectively.  rxx ,  rxv ,  xX2 ,... ,  etc  are  the  Reynolds  stress  tensors. 

In  what  follows,  we  will  reduce  Eq  (6.1)  to  3-D  and  2-D  Euler  equations  for  ihe  convenience  of 
elucidating  the  frequency-domain  POD/ROM  approach. 

6.3.2  Harmonic  Balance 

In  Phase  I,  the  importance  of  aerodynamic  nonlinearities  on  LCO  was  assessed  using  a  novel 
two-dimensional  harmonic  balance  (H.B.)  technique.  With  this  approach,  the  unsteady  flow 
variables  can  be  represented  by  a  Fourier  Series  in  time  with  spatially  varying  coefficients.  This 
assumption  leads  to  a  set  of  harmonic  balance  Euler  equations,  which  can  be  solved  efficiently 
using  conventional  CFD  methods  including  time  marching  with  local  time  stepping  and  multi¬ 
grid  acceleration.  The  two-dimensional  Euler  equations  are: 


au  dF  dG 

dt  dx  dy 


=  0 


where  the  flux  vectors  U,  F  and  G  are  given  by 


(6.2) 
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Next,  representing  conservation  prime  variables  by  sum  of  harmonics  of  fundamental  frequency 
gives: 


P(x,y,t)  =  fu  (x,  y,  t)  -  £  V.  (x,  y), e‘“ 

n  n 

(6A) 

pt  (x,  y,  0  =  X  K  (x’  y)  e'“*  -  Pe  (x>  T,  0  =  Z  (*>  t) 


#/  H 

Requiring  each  frequency  component  to  vanish  independently  (harmonic  balance)  and  collecting 
terms  of  like  harmonics  results  in 


aU(V)  +  3F(V)  +  SOW  .  0  (65) 

dt  dx  8y 

where 

’  U0,  Vq,  E0,  R+ 1,  U+ ,,  V+l,  £+]■■■}  and 

^  =  ico  {•••  0-R0,  0 -U0,  0-V0,  0 -E0,  +1R+1,  +1  -U+l,  +1-F+1,  +1  -£+1  }T  (6.6) 


By  adding  a  ‘pseudo-time’  term,  Eq  (6.5)  can  be  solved  by  a  conventional  CFD  solution,  i.e., 


£V  dF(V)  gG(V)  aU(V) 

dr  dx  dy  dt 


(6.7) 


Note  that  harmonic  balance  equations  are  similar  in  form  to  original  Euler  equations.  Thus, 
existing  Euler  codes  can  be  modified  to  solve  H.B.  equations.  A  two-step  Lax-Wendroff  scheme 
will  be  used.  Also,  since  only  “steady-state”  solution  is  desired,  one  can  use  local  time  stepping, 
multiple-grid  acceleration  techniques  and  residual  smoothing  to  speed  convergence.  The  method 
is  computationally  very  efficient,  at  least  one  to  two  orders-of-magnitude  faster  than  nonlinear 
time-domain  CFD  simulations. 


In  Phase  II,  we  will  use  the  two-dimensional,  nonlinear,  harmonic  balance  aerodynamic  code  and 
add  viscosity  terms,  i.e.,  we  will  develop  a  Navier-Stokes  version  of  this  analysis  suitable  for  the 
analysis  of  oscillating  airfoils.  In  addition,  we  will  use  this  code  to  perform  LCO  analysis  and 
prediction  for  simple  structural  models,  e.g.  plunge  and  pitch  of  a  typical  section  for 
conventional  and  supercritical  airfoils.  In  particular,  we  will  perform  correlation  studies  with  the 
DLR  experimental  data  of  the  NLR  7301  LCO  case. 
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In  Phase  I,  a  three-dimensional  time-linearized  unsteady  aerodynamic  analysis  was  developed. 
In  Phase  II,  this  code  will  be  converted  to  a  harmonic  balance  analysis  of  three-dimensional 
inviscid  transonic  flows.  We  will  use  this  three-dimensional  analysis  to  perform  LCO  analyses 
of  simple  linear  wing  structural  models.  We  will  perform  LCO  correlation  studies  with  a 
business  jet  wing  or  other  wing  to  be  determined. 

In  passing,  we  note  the  following  when  using  aeroelastic  eigenmodes  for  LCO  prediction.  The 
harmonic  balance  method  has  proven  to  be  a  very  effective  method  for  determining  limit  cycle 
oscillations  (LCO)  and  other  nonlinear  responses  of  aeroelastic  systems.  It  is  often  more 
computationally  efficient  and  gives  greater  physical  insight  than,  for  example,  time  marching 
simulations.  On  the  other  hand  it  has  some  limitations,  e.g.  it  does  not  allow  an  investigation  of 
the  possible  dependence  of  the  response  on  the  initial  conditions.  Nor  does  it  allow  a  precise 
assessment  of  non-periodic  and  chaotic  motions.  At  best  the  harmonic  balance  method  can  only 
approximate  a  non-periodic  motion.  In  practice  these  limitations  are  not  often  serious 
drawbacks;  however  these  observations  do  suggest  that  future  work  may  well  involve  a  balanced 
consideration  of  both  time  simulations  and  frequency  domain  (harmonic  balance)  studies. 

Another  consideration  is  the  number  of  harmonics  to  be  retained  in  a  harmonic  balance  analysis. 
If  there  is  a  dominant  harmonic,  a  single  harmonic  may  suffice,  although  it  will  be  important  to 
assess  the  effects  of  higher  harmonics  on  the  fundamental  harmonic  per  se.  And  finally  the 
algebraic  complexity  of  using  the  method  for  multi-modal  systems  may  quickly  get  out  of  hand 
unless  special  measures  are  taken  to  deal  with  this  issue.  We  will  return  to  this  latter  issue  as, 
even  with  reduced  order  aerodynamic  modeling  of  the  structure  and  fluid,  the  order  of  the 
reduced  model  will  be  larger  than  that  often  encountered  in  classical  harmonic  balance  analyses, 
i.e.,  on  the  order  of  several  tens  of  states. 

6.3.3  Time-Linearized  Flow  Three-Dimensional  Flow  Solver 


In  Phase  I,  we  developed  a  three-dimensional  time-linearized  Euler  analysis  of  unsteady  flow 
abo'\  a  wing.  We  start  with  the  assumption  that  the  unsteady  flow  u (x,y,z,t)  may  be  expanded 
in  a  perturbation  series  of  the  form 


u(x,  y,  z,  t)  =  U(x,  y,  z,  t )  +  u(x,  y,  z)e,cot 


(6.8) 


where  U {x,y,z,t)  represents  the  steady  background  flow,  and  u(x,y,z)  is  the  complex  amplitude 
of  the  unsteady  small-disturbance  flow  arising  from  the  wing  vibration,  which  vibrates  at 
frequency  co.  Substituting  Eq  (6.8)  into  the  nonlinear  three-dimensional  Euler  equations,  and 
expanding  the  result  in  a  perturbation  series  in  the  small-disturbance  quantities,  one  finds  that  to 
zeroth  and  first  order  the  governing  equations  reads  respectively 


^(U)  <3G(U)  £1(U) 

dx  dy  dz 


(6.9) 
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(6.10) 


d(  cF  \  d(dG  \  dfcH. 

lC0U+  +  4>{&juJ  +  &Va]  UJ 

In  a  conventional  time-linearized  analysis,  Eq  (6.9)  and  Eq  (6.10)  are  solved  using  pseudo  time 
marching  and  standard  CFD  techniques.  Note  that  since  time  does  not  appear  explicitly  in  either 
equation,  they  both  may  be  solved  using  computational  acceleration  techniques  such  as  multiple 
grid  acceleration,  local  time  stepping,  and  residual  smoothing.  The  result  is  that  these  two 
equations  can  be  solved  at  a  fraction  of  the  cost  of  solving  the  full  nonlinear  Euler  equations 
using  conventional  time  marching. 

In  Phase  I,  a  simple  clean  wing  model  was  developed.  This  CFD  analysis  used  a  simple  grid 
structure  that  limited  the  approach  to  clean  wings.  In  Phase  II,  this  model  will  be  extended  to 
allow  more  complicated  grid  topologies,  which  will  allow  us  to  analyze  multi-body 
configurations,  e.g.  wing/store. 

6.3.4  Reduced  Order  Modeling  of  Time-Linearized  Aerodynamic  Models 

Having  developed  a  three-dimensional  time-linearized  flow  solver,  we  next  consider  the 
reduction  of  that  model  using  proper  orthogonal  decomposition  techniques  (POD).  The  idea 
behind  the  frequency  domain  proper  orthogonal  decomposition  is  a  simple  one.  We  first 
calculate  the  small-disturbance  response  of  the  aerodynamic  system  at  M  different  combinations 
of  frequency  and  excitation.  The  solutions,  also  called  “snapshots,”  are  denoted  by  qm  for 
m-1,2, ... ,M .  These  snapshots  are  then  linearly  combined  to  for  a  smaller  number  of  basis  vectors 
<f>k  for  k=l,2,  ...,K  where  K<M.  In  other  words, 

(j>k  =Qvk  for  k  =  1,2,...,  K  (6.11) 

where  the  mlh  column  of  Q  is  just  qm.  In  the  proper  orthogonal  decomposition  technique,  the 
vectors  vk  are  found  by  solving  a  small  eigenvalue  problem  of  the  form 

Q"Qv,=V*  (6A2) 

Only  the  eigenvectors  vk  v/iih  the  largest  eigenvalues  Xk  are  used  to  form  basis  vectors  defined 
by  Eq  (6. 1 1).  Qw  is  the  Hermitian  of  Q. 

Having  computed  the  POD  basis  vectors,  we  assume  that  they  will  provide  a  useful  basis  for 
computing  the  unsteady  solution  at  some  other  frequency  and/or  external  excitation  than  was  use 
to  general  the  original  snapshots.  Thus,  we  let 

q  =  <t»£  (6.13) 

where  O  is  an  NxK  matrix  whose  kth  column  is  simply  the  basis  vector  4>k,  and  ^  is  a  vector  of 
aerodynamics  state  variables. 

We  note  that  when  discretized,  Eq  (6.10)  has  the  form 


Aq  =  A0q  +  z<yA,q  =  b0  +  ico b, 


(6.14) 
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where  A0  and  A|  are  independent  of  the  excitation  frequency  co,  and  are  purely  real,  and  q  is  the 
vector  containing  u  at  all  the  nodes  of  the  CFD  grid. 

Substitution  of  Eq  (6. 1 1)  into  Eq  (6.14),  and  projection  of  the  error  onto  the  space  spanned  by  the 
basis  vectors  gives 

OwAO£=  Ayj£=  Owb  (6.15) 

Finally,  the  reduced-order  aerodynamic  matrix  Ar,  which  is  much  smaller  than  the  original 
aerodynamic  matrix  A,  is  factored  using  LU  decomposition,  and  Eq  (6.15)  is  solved  for  the 
unknown  aerodynamic  state  variables.  This  step  is  computationally  very  efficient  because  the 
reduced-order  aerodynamic  matrix  is  small,  sometimes  as  small  as  10x10,  but  rarely  larger  than 
100x100.  The  major  expense  in  constructing  the  reduced-order  aerodynamic  model  is  the 
computation  of  the  snapshots;  the  cost  of  finding  the  basis  vectors  and  solution  to  Eq  (6.15)  is 
negligible  by  comparison. 

Having  described  the  basic  reduced-order  modeling  technique,  we  next  describe  how  to 
incorporate  an  aerodynamic  reduced-order  model  into  an  aeroelastic  model  of  flutter.  Consider, 
for  example,  a  two  degree-of-ffeedom  structural  dynamic  model  of  a  typical  section.  The 
governing  equations  of  motion  are  of  the  form 

Mh  +  Kh  =  f  where  h  =  { h,  a}T  (6. 1 6) 

and  h  and  a  are  the  plunging  and  pitching  degrees  of  freedom  of  the  typical  section.  M  and  K  are 
the  mass  and  stiffness  matrices,  respectively. 

Note  that  the  aerodynamic  force  vector  f  is  obtained  from  integrals  involving  the  pressure  at  the 
surface  of  the  airfoil.  When  discretized,  these  integrals  may  be  expressed  as 

f  =  C  q  (6.17) 

where  C  is  a  sparse  2xN  matrix.  Similarly,  for  the  case  of  airfoil  vibration,  the  vector  b  on  the 
right-hand  side  of  Eq  (6.14)  can  be  expressed  as 


b  =  b0  +  iabl  =  B0h  +  i  «B,  h 


(6.18) 


where  now  we  have  made  the  assumption  that  the  airfoil  motion  is  harmonic  in  time  although  co 
may  be  complex).  For  large  CFD  models,  finding  the  eigenvalues  is  prohibitively  expensive.  To 
reduce  the  size  of  the  model,  we  again  assume  that  the  number  of  aerodynamic  states  can  be 
reduced  using  Eq  (6. 1 3).  Again,  projecting  the  error  of  the  aerodynamic  equations  onto  the  space 
spanned  by  the. aeroelastic  basis  vectors  gives  the  desired  reduced-order  aeroelastic  model,  i.e., 
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(6.19) 
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6.3.5  Reduced-Order  Aeroelastic  Model  for  Multi-Degree-Of-Freedom  Structures 

In  Eq  (6.19),  we  have  derived  a  reduced-order  aeroelastic  model  for  a  two  degree  of  freedom  (2 

d.o.f.)  structure,  i.e.,  h  =  {h,a}T .  For  multi-degree-of-freedom  structures,  h  can  also  be 
interpreted  as  the  displacements  of  each  degree  of  freedom  in  the  structural  model.  In  this 
situation,  the  size  of  Eq  (6.19)  becomes  K+2n,  where  n  is  the  number  of  structural  degrees  of 
freedom.  For  a  realistic  wing  structural  model,  n  can  be  on  the  order  of  thousands,  rendering  Eq 
(6.19)  a  very  large  size  eigenvalue  problem.  Solving  such  a  large  eigenvalue  problem  would 
practically  be  impossible. 

The  Modal  Approach 

To  circumvent  this  problem,  we  introduce  the  so-called  “ modal  approach ”  to  Eq  (6.19).  The 
modal  approach  approximates  h  by  the  superposition  of  the  low-order  structural  modes,  i.e., 

h  =  '¥arj  (6-2°) 

where  To  is  the  modal  matrix  whose  columns  contain  the  modal  data  of  the  low-order  structural 
modes  and  r\  are  is  the  generalized  coordinate  vector.  Since  the  magnitude  of  the  modes  can  be 
arbitrary,  they  are  usually  normalized  by  the  square  root  of  their  respective  generalized  masses, 
giving  a  unit  generalized  mass  matrix.  The  justification  for  using  the  modal  approach  is  based 
on  the  premise  that  most  of  the  aeroelastic  responses  are  dominated  by  the  lower-order  structural 
modes.  Usually,  for  a  single  wing  structure,  the  lowest  ten  (10)  structural  modes  are  sufficient  to 
accurately  represent  h.  Substituting  Eq  (6.20)  into  Eq  (6.19)  and  pre-multiplying  the  second  and 

third  rows  of  Eq  (6. 1 9)  by  Tj  yields 


ARo  -OX*, 
0  0 
-VJC0  [»,2] 


>  +  10) 
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f*l 

I77 
l  }l 


=  0 


(6.21) 


where  co„  and  r^  are  the  natural  frequency  and  the  generalized  mass  of  the  zth  structural  mode, 

respectively.  Now,  the  size  of  Eq  (6.21)  becomes  K+2m,  where  m  is  the  number  of  structural 
modes.  Because  m  is  much  less  than  the  number  of  structural  d.o.f.  »,  the  size  of  the  aeroelastic 
system  is  substantially  reduced. 


The  Structural/Modal  Damping 

Note  that  a  structural/modal  damping  matrix 


has  been  added  to  Eq  (6.21). 


With  this 


added  matrix,  Eq  (6.21)  will  facilitate  our  subsequent  study  of  LCO.  In  an  earlier  LCO  study, 
ZONA  has  suggested  that  the  structural  damping  could  play  a  decisive  role  in  LCO  for  a 
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wing/store  system.  For  LCO  study  using  the  proposed  method,  one  could  alternatively  use  Eq 
(6.21)  to  delineate  the  effects  with  and  without  structural  damping.  Also  note  that  the  modal 
approach  in  fact  increases  the  sparseness  of  the  matrices  in  Eq  (6.21).  Thus,  the  eigenvalue 
problem  of  Eq  (6.21)  can  be  solved  more  efficiently  by  a  sparse  eigensolver. 

6.3.6  Structural  Compatible  Reduced-Order  Aeroelastic  Model 

One  of  the  technical  issues  involved  in  the  CFD  aeroelastic  computations  is  the  displacement 
transferal  from  the  structural  finite  element  grid  to  the  CFD  surface  grid.  This  technical  issue 
arises  from  the  problem  where  the  CFD  surface  grid  and  the  structural  finite  element  (FEM)  grid 
are  considerably  different;  in  their  locations  and  densities.  Solving  such  a  displacement 
transferal  problem  of  complex  configuration  such  as  whole  aircraft  with  external  stores  is  by  no 
means  a  trivial  task. 

In  March  1999,  ZONA  received  a  two-year  contract  from  NASA/Langley  (Ref  21)  to  develop  a 
Boundary  Element  Method  (BEM),  called  the  BEM  Solver,  for  the  data  transferal  between  the 
FEM  grid  and  the  CFD  surface  grid  (see  Section  6.1,  ZONA’s  Related  Work).  By  formulating 
the  data  transferal  problem  as  an  equivalent  solid  mechanics  problem,  the  BEM  Solver  generates 
a  universal  spline  matrix  [S]  which  relates  the  displacements  at  the  FEM  grid  to  the  CFD  grid 
such  that 

=  S  ¥,  (6.22) 

where  and  are,  respectively,  the  modal  matrix  at  the  CFD  grid  and  at  the  FEM  grid. 
With  the  universal  spline  matrix  [S]  at  hand,  converting  Eq  (6.21)  to  a  structural-compatible 
reduced-order  aeroelastic  model  is  straightforward.  Substituting  Eq  (6.22)  into  Eq  (6.21)  gives: 

aRo  -d>"B0s^ 


Lm<  J 

_-v,TsTc«>  k,!]  k»,<r] 

Finally,  we  arrive  at  the  structural-compatible,  modal-based,  reduced-order  aeroelastic  model,  Eq 
(6.23).  The  size  of  this  model  is  K+2m,  where  m  is  the  number  of  low-order  structural  modes 
and  is  usually  on  the  order  of  ten  for  a  single  wing  structure.  Eq  (6.23)  contains  two  very  sparse 
matrices  that  can  be  solved  efficiently  by  a  sparse  eigenvalue  solver. 
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=  0  (6.23) 


6.3.7  Modeling  LCO  of  High  Degree-of-Freedom  Nonlinear  Systems 

The  nonlinear  system  that  we  wish  to  model  can  for  the  most  part  be  modeled  as  quasi-linear 
nonlinear  vector  equations  of  the  form 
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ztyMq-t-N(q)=  0 


(6.24) 


where  q  is  a  very  large  matrix  arising  containing  the  Fourier  coefficients  of  the  unknowns  in  a 
harmonic  balance  analysis,  M  is  the  a  constant  matrix,  and  N  is  a  nonlinear  vector  operator 
arising  from  the  harmonic  balance  analysis.  Here  q  would  contains  the  unknown  flow  solution  in 
the  three-dimensional  harmonic  balance  of  the  three-dimensional  flow  field,  and  also  the 
structural  dynamic  modes.  Thus,  the  system  may  contain  hundreds  of  thousands  of  degrees  of 
freedom.  This  system  of  equations  is  solved  using  pseudo-time  time  marching.  Thus,  Eq  (6.24) 
is  solved  by  marching  the  equation 


—  +  ityMq  +  N(q)=  0  (6-25) 

dt 

in  time  until  a  steady  state  is  reached.  However,  if  co  is  not  known  accurately,  then  Eq  (6.25) 
will  not  converge,  but  will  itself  go  into  a  mathematical  limit  cycle.  One  can  show,  however, 
that  the  solution  q  will  be  nearly  correct.  The  solution  can  be  improved  by  first  computed  a 
better  estimate  for  musing  a  nonlinear  Rayleigh  quotient,  i.e. 


q"N(q) 

qwM"M(q) 


(6.26) 


Eq  (6.25)  can  then  be  marched  again  to  improve  the  estimate  of  q.  This  process  can  be  repeated 
until  convergence.  The  result  is  a  description  of  the  LCO  behavior  of  a  very  high  dimensional 
system,  i.e.  a  nonlinear  CFD  model  coupled  to  a  linear  or  nonlinear  structural  model.  If  Phase  II, 
this  technique  will  be  applied  to  the  harmonic  balance  flow  solver  coupled  to  a  linear  and/or 
nonlinear  structure. 


6.3.8  Selection  of  Test  Cases 

Test  cases  for  the  validation  of  the  proposed  methodology  include: 

•  NLR  7301  supercritical  airfoil  transonic  LCO  case 

•  MAVRIC-I  (Aeroelastic  Validation  Research  Involving  Computations)  business  jet  wing 
LCO  case 

•  Three  F- 1 6/Store  LCO  cases 

•  Five  F/A- 1 8/Store  LCO  cases 

NLR  730 J  Supercritical  Airfoil  Transonic  LCO  Case 

The  details  of  this  case  has  already  been  discussed  in  Chapter  5.  This  case  is  selected  to  validate 
the  2-D  nonlinear  harmonic  balance  aerodynamic  code  using  the  Euler  equations  and  the  Navier- 
Stokes  equations.  The  accuracy  of  the  2-D  nonlinear  harmonic  balance  aerodynamic  code  will 
be  assessed  by  comparing  its  predicted  LCO  frequency  and  amplitude  with  the  DLR 
experimental  data  whereas  the  efficiency  will  be  assessed  by  comparing  its  computer  time  with 
the  CFL3D  time-marching  computations. 
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MAVRIC-I  Business  Jet  Wing  LCO  Case 


MAVRIC-I  business  jet  wing  is  an 
ongoing  wind  tunnel  test  program  of 
NASA/Langley  aeroelasticity  branch. 

The  model  is  constructed  with  a  simple 
stepped  aluminum  plate  providing  the 
wing  stiffness  and  fitted  with  end-grain 
balsa  wood  to  provide  the  wing 
contour.  Previous  experiment  has 
shown  that  in  the  0.8  to  0.9  Mach 
number  range,  the  model  motions  were 
predominantly  in  the  first  wing-bending 
mode  and  exhibited  LCO. 

Currently,  the  model  is  being  re-tested  in  the  Transonic  Dynamic  Tunnel  (TDT)  as  the 
MAVRIC-I.  Fig  6.3  indicates  the  location  of  instrumentation  that  has  been  added  to  the  model. 
Three  chords  of  unsteady  pressure  transducers  are  installed  at  span  stations  0.22,  0.63  and  0.87. 
Each  chord  has  28  upper  surface  and  18  lower  surface  close-mounted  transducers.  Eight 
accelerometers  are  mounted  along  the  leading  and  trailing  edge  of  the  wing  and  bending/torsion 
strain  gages  are  installed  at  the  root.  The  intent  of  the  retest  is  to  obtain  unsteady  pressure  and 
wing  response  data  under  conditions  of  LCO  in  order  to  validate  CFD  codes  for  such  conditions. 
This  retest  program  will  be  complete  within  6-9  months. 

Because  of  its  simple  finite  element  modeling  and  simple  aerodynamic  geometry,  the  MAVRIC- 
I  business  jet  wing  model  is  an  ideal  LCO  test  case  for  validating  the  3-D  harmonic  balance  code 
for  single  body  (to  be  developed  in  the  first  year  of  Phase  II).  Because  all  harmonics  of  the 
measured  time-domain  unsteady  pressures  at  LCO  can  be  easily  obtained  by  a  Fourier  series 
techniques,  a  detailed  pressure  distribution  comparison  between  the  measurements  and  the 
prediction  by  3-D  harmonic  balance  code  in  terms  of  frequency-domain  harmonics  can  be 
performed.  Of  course,  the  comparison  of  LCO  frequency  and  amplitude  at  various  Mach 
numbers  and  dynamic  pressures  will  also  be  conducted. 

Three  F-l  6/Store  LCO  Cases 

Through  a  collaboration  project  with  Lockheed  Martin  Tactical  Aircraft  System  and  Eglin  Air 
Force  Base  (see  Section  6.1,  ZONA’s  Related  Work),  ZONA  has  established  a  large  database  of 
the  F-l  6/store  configurations  including  their  finite  element  models  and  flight  test  LCO  data.  In 
general,  the  F-l 6  aircraft  with  various  stores  has  experienced  three  types  of  LCO;  (1)  LCO  starts 
from  M  =  0.6  and  disappears  at  M  >  1.0.  This  case  is  denoted  as  “ F-l  6  LCO  Type  A ”.  (2)  LCO 
appears  only  in  the  transonic  flight  regime.  This  case  is  denoted  as  “F-l  6  LCO  Type  B". 
(3)  LCO  starts  from  the  transonic  Mach  numbers  and  sustains  into  the  supersonic  flight  regime. 
This  case  is  denoted  as  “F-l 6  LCO  Type  C\ 


Pressure  Trans  Aicer  Rows 


Figure  6.3  Instrumentation  Layout  for 
Refurbished  MAVRIC-I  Business  Jet  Wing  Model 
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Figure  6.4  Three  Selected  F-16/Store  LCO  Cases 

In  order  to  validate  the  3-D  harmonic  balance  code  with  a  multi-body  capability  (to  be  develop  in 
the  second  year  of  Phase  II),  three  F- 16/store  LCO  cases  are  selected;  each  represents  one  type 
of  LCO  discussed  above.  These  three  F- 16/store  configurations  and  their  respective  LCO  flight 
conditions  are  shown  in  Fig  6.4. 


Five  Selected  F/A-18/Store  LCO  Cases 

In  1999,  ZONA  has  received  a  Navy  SBIR  Phase  I  contract  entitled  “. Adaptive  Reconfigurable 
Control  Based  on  a  Reduced-Order  System  Identification  for  Flutter  and  Aeroservoelastic 
Instability  Suppression"  (see  Section  6.1,  ZONA’s  Related  Work).  As  one  of  the  team  members, 
the  Boeing  company,  St.  Louis,  supplied  ZONA  with  all  the  necessary  data  of  the  F/A-18  LCO 
configurations  including  the  structural  finite  element  models  with  various  stores,  the  LCO  flight 
test  data  and  flight  control  system  for  F/A-18  LCO  suppression.  In  general,  two  types  of  LCO 
frequencies  were  observed  during  the  F/A- 18/store  flight  tests:  (1)  Type  A:  LCO  frequency  at 
5.6  Hz  for  cases  of  wing/store  with  tip  missile,  (2)  Type  B:  LCO  frequency  at  8.8  Hz  for  cases 
of  wing/store  with  tip  launcher  only. 
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Case 

Store  Configuration 

LCO  Frequency 

Flight  Condition 

1. 

- * 

Wing  Tip:  Launcher  +  Missile 
Outboard  Pylon:  MK-84 

Inboard  Pylon:  None 

5.6  Hz 

M  *  0.88  -  1.0 

2. 

Wing  Tip:  Launcher  +  Missile 
Outboard  Pylon:  MK-84 

Inboard  Pylon:  MK-84 

5.6  Hz 

M  «  0.88-  1.0 

#  • 

3. 

Wing  Tip:  Launcher  +  Missile 
Outboard  Pylon:  AGM-88 
Inboard  Pylon:  MK-84 

5.6  Hz 

M  «  0.88-  1.0 

j  *  *  * 

4. 

Wing  Tip:  Launcher  Only 
Outboard  Pylon:  MK-83 

Inboard  Pylon:  None 

8.8  Hz 

M  >  0.9 

* 

5. 

Wing  Tip:  Launcher  Only 
Outboard  Pylon:  MK-83 

Inboard  Pylon:  MK-83 

8.8  Hz 

M  >  0.9 

» 

i 

i 

Figure  6.5  Five  Selected  F/A-18  /  Store  LCO  Cases 


In  order  to  further  validate  the  3-D  harmonic  balance  code  with  a  multi-body  capability,  we 
selected  five  F/A- 18/store  LCO  test  cases  (three  for  type  A  and  two  for  type  B)  whose 
configuration  and  flight  conditions  are  presented  in  Fig  6.5. 

To  finish  the  three  F- 16/store  cases  and  five  F/A- 18/store  cases  within  a  two-year  project  seems 
to  be  a  very  ambitious  plan.  However,  this  is  not  the  case  if  the  POD/ROM  methodology  is 
employed.  Since  it  is  known  that  the  underwing  stores  generally  has  little  aerodynamic  influence 
on  the  aeroelastic  characteristics,  i.e.,  only  the  inertial  effects  of  the  underwing  stores  impact  the 
aeroelastic  characteristics,  aeroelastic  analysis  of  wing/store  configurations  can  exclude  the 
underwing  stores  in  the  aerodynamic  modeling  but,  of  course,  include  their  inertial  effects  in  the 
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structural  modeling.  Based  on  this  assumption,  the  three  F-16/store  and  five  F/A-l 8/store  cases 
become  the  ideal  test  cases  for  the  POD/ROM  methodology.  As  already  demonstrated  in  our 
Phase  I  achievements  presented  in  Chapter  4,  one  can  obtain  accurate  POD/ROM  unsteady 
aerodynamics  using  solution  snapshots  unrelated  to  the  actual  structural  modes.  This  suggests 
that  only  four  sets  of  POD/ROM  computations  are  required,  two  for  F-16  with  and  without  tip 
missiles  and  another  two  for  F/A-l 8  with  and  without  tip  missiles,  but  using  simple  wing  motion 
snapshots  that  are  unrelated  to  the  actual  F-16  or  F/A-l 8  structural  modes.  Once  these  snapshot 
solutions  are  at  hand,  the  unsteady  aerodynamics  of  the  three  F-l  6/store  and  five  F/A-l  8/store 
cases  can  be  obtained  by  using  the  snapshot  ensemble  technique  described  in  Chapter  4,  but  with 
respect  to  their  actual  structural  modes  and  tip  missile/launcher  aerodynamic  configurations. 

Note  that  excluding  the  underwing  stores  from  the  aerodynamic  modeling  is  based  on  the 
assumption  of  their  little  aerodynamic  influence  on  the  whole  aircraft  aeroelastic  characteristics. 
To  verify  this,  we  will  conduct  a  flutter  analysis  on  one  of  the  wing/store  configurations  but 
including  all  underwing  stores  in  the  aerodynamic  modeling.  The  flutter  result  will  be  used  to 
verify  the  snapshot  procedure  discussed  above.  Once  verified,  this  implies  that  the  POD/ROM 
methodology  can  be  developed  as  an  highly  efficient  computational  tool  for  massive  flutter  and 
LCO  analyses  of  aircraft  with  various  underwing  stores. 

6.3.9  Statement  of  Work 

Tasks  to  be  performed  in  Phase  II  are  defined  below.  The  schedule  for  the  Phase  II  statement  of 
work  is  shown  in  Section  6.7.10.  Timing  for  related  program  deliverables,  meeting  and 
presentations  are  also  noted  on  the  schedule. 

Task  1:  Development  of  a  2-D  Navier-Stokes  Solver  with  Harmonic  Balance  (NSHB)  Code 

Take  the  two-dimensional,  nonlinear,  harmonic  balance  aerodynamic  code  using  the  Euler 
equations  which  was  developed  in  Phase  I  and  add  viscosity,  i.e.,  this  becomes  a  Navier-Stokes 
code. 

Task  2:  Validation  of  the  2-D  NSHB  Code  for  Simple  Structures  in  Transonic  Flow 

Use  the  2-D  NSHB  code  to  perform  LCO  analysis  and  prediction  for  simple  structural  models, 
e.g.,  plunge  and  pitch  of  typical  sections  for  conventional  and  supercritical  airfoils  in  transonic 
flow. 

-  Perform  flutter  and  LCO  studies  of  the  NACA  640 10A  conventional  airfoil  described  in 
Chapter  3.  Assess  the  impact  of  viscosity  on  transonic  flutter  dip  and  LCO  by  comparing 
the  results  of  the  2-D  NSHB  with  those  of  the  inviscid  2-D  harmonic  balance  code 
developed  in  Phase  I. 

-  Perform  flutter  and  LCO  studies  of  the  NLR  7301  supercritical  airfoil  described  in 
Chapter  5.  Correlate  the  predicted  LCO  frequency  and  amplitude  with  the  DLR 
experimental  results  and  assess  the  efficiency  of  the  2-D  NSHB  code  by  comparing  its 
computer  time  with  that  of  the  time-marching  CFL3D  computation. 

Task  3:  Development  of  the  3-D  Euler  Solver  with  the  Harmonic  Balance  (EHB)  Code 
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Take  the  3-D  time-linearized  aerodynamic  code  using  the  Euler  equations  developed  in  Phase  I 
and  convert  it  to  a  fully  nonlinear  dynamics  code  using  the  harmonic  balance  method  previously 
developed  for  the  2-D  code. 

Task  4:  Validation  of  the  3-D  EHB  Code  for  Simple  Wing  Structural  Model 

Use  the  3-D  EHB  code  to  perform  LCO  analysis  and  prediction  for  simple  wing  structural 
models.  Perform  correlation  studies  with  the  MAVRIC-I  business  jet  wing. 

-  Generate  a  3-D  mesh  for  the  MAVRIC-I  business  jet  wing. 

-  Use  the  BEM  Solver  to  transfer  the  mode  shapes  computed  at  the  structural  finite  element 
grid  to  the  CFD  surface  mesh  of  the  MAVRIC-I  business  jet  wing. 

-  Perform  LCO  analysis  at  selected  transonic  Mach  numbers  and  correlate  the  predicted 
frequency-domain  pressure  distribution  with  the  NASA/Langley  TDT  measurements.  Assess 
the  accuracy  of  the  predicted  LCO  frequency  and  amplitude  by  the  comparison  with  the  TDT 
measurements. 

Task  5:  Development  of  a  2-D  NSHB  Code  with  Multi-Structural  Mode  Capability 

Take  the  2-D  NSHB  code  developed  in  Task  1  and  add  a  multi-structural  mode  capability.  This 
will  provide  the  capability  for  LCO  analysis  of  complex  multi-mode,  nonlinear  structural  models 
such  as  an  airfoil  and  control  surface  with  ffeeplay. 

Task  6:  Validation  of  the  2-D  NSHB  Code  for  Nonlinear  Aerodynamic  and  Structural 
Models 

Convert  the  2  degrees  of  freedom  NACA  640 10A  and  NLR  7301  airfoils  to  a  3  degrees  of 
freedom  structure  by  adding  a  control  surface  motion  with  ffeeplay.  Use  the  2-D  NSHB  code 
with  multi-structural  mode  capability  to  perform  LCO  studies  of  the  3  degrees  of  freedom 
NACA  640 10A  and  NLR  7301  airfoils.  Establish  a  test  bed  for  novel  and  improved  harmonic 

balarce  solution  methods  to  predict  LCO  with  nonlinear  aerodynamics  and  structural  models. 

> 

Task  7:  Development  of  a  3-D  Time-Linearized  Euler  Solver  (TLE)  Code  with  a  Multi- 
Body  Capability 

Take  the  3-D  time-linearized  aerodynamic  code  using  the  Euler  equations  developed  in  Phase  I 
and  add  a  multi-body  capability.  This  will  provide  the  capability  to  predict  the  transonic  flutter 
dip  of  wing  with  stores  configurations. 

Task  8:  Validation  of  the  3-D  TLE  Code  for  Multi-Body  Configurations 

Use  the  3-D  TLE  code  with  the  multi-body  capability  to  perform  flutter  analysis  of  the  three  F- 
1 6/store  and  five  F/A-  18/store  cases.  Correlate  the  predicted  flutter  boundaries  with  the  flutter 
flight  test  data. 

-  Generate  the  multi-block  mesh  systems  of  the  F- 16/store  and  F/A- 18/store  configurations 
with/without  the  tip  missile. 
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-  Use  the  POD/ROM  methodology  to  obtain  the  snapshot  solutions  of  the  F- 16/store  and  F/A- 
1 8/store  configurations  with  simple  wing  motions  that  are  unrelated  to  the  actual  structural 
modes. 

-  Transfer  the  structural  mode  shapes  of  the  three  F- 16/store  and  five  F/A- 18/store  cases  from 
the  structural  finite  element  grid  to  the  CFD  surface  grid  using  the  BEM  solver. 

-  Compute  the  unsteady  aerodynamics  of  these  F- 16/store  and  F/A- 18/store  cases  with  respect 
to  their  actual  structural  modes  using  the  snapshot  ensemble  technique. 

-  Perform  flutter  analysis  of  the  three  F-  16/store  and  five  F/A- 18/store  cases.  Correlate  the 
predicted  flutter  boundaries  with  the  flutter  flight  test  data. 

Task  9:  Development  of  a  3-D  Harmonic  Balance  Code  with  Multi-Body  Capability 

Take  the  3-D  TLE  code  developed  in  Task  7  and  add  the  harmonic  balance  method.  This  will 

provide  the  capability  of  LCO  analysis  for  multi-body  geometry  like  the  wing/store 

configurations. 

Task  10:  Development  of  a  Rayleigh  Quotient  Technique  for  LCO  Prediction  of  High 
Degree-of-Freedom  Nonlinear  System 

Apply  Rayleigh  Quotient  technique  to  the  codes  developed  in  Tasks  2,  3,  or  9  to  demonstrate 
feasibility  of  computing  LCO  with  high  degree  of  freedom  nonlinear  system.  This  will  give  a 
description  of  the  LCO  behavior  of  a  very  high  dimensional  system,  i.e.,  a  nonlinear  CFD  model 
coupled  with  a  nonlinear  structural  model. 

Task  11:  Validation  of  the  3-D  Harmonic  Balance  Code  for  the  Multi-Body  Configurations 

Use  the  3-D  harmonic  balance  code  with  the  multi-body  capability  developed  in  Task  9  to 

perform  LCO  analysis  and  prediction  of  the  three  F- 16/store  and  five  F/A- 18/store  cases. 

Correlate  the  predicted  LCO  frequencies  and  amplitudes  of  these  cases  with  their  respective 
flight  test  data. 
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6.3.10  Planned  Program  Schedule 


Yr  1  Quarter 
1  I  2  I  3  I  4 


Yr  2  Quarter 
1  I  2  |  3  I  4 


Performed  by 


1 .  Development  of  a  2-D  Navier-Stokes  Solver  with . 

Harmonic  Balance  (NSHB)  Code 

2.  Validation  of  the  2-D  NSHB  Code  for  Simple 
Structures  in  Transonic  Flow 

-  Perform  flutter  and  LCO  studies  of  the  NACA . 

640 10A  conventional  airfoil 

-  Perform  flutter  and  LCO  studies  of  the  NLR . 

7301  supercritical  airfoil 

3.  Development  of  the  3-D  Euler  Solver  with  the . 

Harmonic  Balance  (EHB)  Code 

4.  Validation  of  the  3-D  EHB  Code  for  Simple  Wing 
Structural  Model 

-  Generate  a  3-D  mesh  for  the  MAVRIC-I . 

business  jet  wing 

-  Use  the  BEM  Solver  to  transfer  the  mode  shapes 
from  the  structural  grid  to  the  CFD  surface  grid 

-  Perform  LCO  analysis  at  transonic  Mach . 

numbers  of  the  MAVRIC-I  business  jet  wing 

5.  Development  of  a  2-D  NSHB  Code  with  Multi- . 

Structural  Mode  Capability 

6.  Validation  of  the  2-D  NSHB  Code  for  Nonlinear . 

Aerodynamic  and  Structural  Models 

7.  Development  of  a  3-D  Time-Linearized  Euler . 

Solver  (TLE)  Code  with  a  Multi-Body  Capability 

8.  Validation  of  the  3-D  TLE  Code  for  Multi-Body 
Configurations 

-  Generate  the  multi-block  mesh  systems  of  the  F- 
1 6/store  and  F/A- 18/store  configurations 

-  Use  the  POD/ROM  methodology  to  obtain  the . 

snapshot  solutions 

-  Transfer  the  structural  mode  shapes  from  the . 

structural  finite  element  griH  to  the  CFD  surface 
grid  using  the  BEM  solver 

-  Compute  the  unsteady  aerodynamics  of  these  F-  " 
16/store  and  F/A- 18/store  cases  using  the 
snapshot  ensemble  technique 

-  Perform  flutter  analysis  of  the  three  F- 16/store . 

and  five  F/A-  18/store  cases 

9.  Development  of  a  3-D  Harmonic  Balance  Code 
with  Multi-Body  Capability 

10.  Development  of  a  Rayleigh  Quotient  Technique  for 
LCO  Prediction  of  High  Degree-of-Freedom 
Nonlinear  System 

1 1 .  Validation  of  the  3-D  Harmonic  Balance  Code  for . 

the  Multi-Body  Configurations 

12.  Documentation 

-  Interim  Report 

Final  Report . 

Kick-Off  Meeting . 

Final  Oral  Presentation . 


ZONA/Duke 


ZONA 


ZONA 


ZONA 


Duke/ZONA 


ZONA 


ZONA 


Duke/ZONA 


ZONA 


ZONA/Duke 


ZONA/Duke 

ZONA/Duke 


6.3.11  Deliverables 


The  deliverables  and  delivery  cycle  of  Phase  II  are  presented  in  Table  6. 1 . 


Table  6.1  Deliverables  of  the  Proposed  Phase  II  Project 


Deliverable 

Delivery  Cycle 

Interim  Reports 

Quarterly 

Funds  &  Man  Hour  Expenditure  Reports 

Quarterly 

Final  Report 

As  Required 

Executable  Code  of  the  Final  Program 

As  Required 
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CHAPTER  7 


COMMERCIALIZATION  STRATEGY 


With  a  strong  technical  background,  ZONA  also  has  extensive  experience  in  commercialization 
of  its  product.  In  fact,  ZONA  has  been  serving  the  aerospace  community  through 
consulting/contractual  work  and  software  licensing  support  since  1986. 

•  ZONA’s  unsteady/steady  aerodynamic  software  in  the  aerospace  market  most  notably 
includes  ZONA51  (currently  the  Aero  Option  II  in  MSC/NASTRAN  -  an  industry  standard) 
and  the  ZAERO  software  system  (covering  the  entire  Mach  range).  ZONA  codes  have,  thus 
far,  accumulated  over  120  users  worldwide. 

•  Under  an  AF/STTR  Phase  II  contract,  ZONA  and  MSC  Software  are  reaching  a  business 
agreement  to  commercialize  and  jointly  market  ASTROS*  (the  seamlessly  integrated 
ZAERO  module  into  ASTROS).  ASTROS  is  a  popular  Automated  STRuctural  Optimization 
System  software  among  aerospace  industry  and  universities. 

7.1  Next  Generation  Aeroelastic  Software 

After  some  20  years  of  continuing  R&D  effort  in  unsteady  transonic  aerodynamics,  high-level 
(3-D  Euler  and  Navier-Stokes)  unsteady  CFD  methods  remain  inadequate  to  handle  the  fixed- 
wing  transonic  aeroelastic  problems  encountered  in  industry.  ZONA’s  recent  marketing  survey 
reveals  the  possible  causes  of  such  an  inadequacy.  These  include  that:  (1)  the  current  CFD 
methods  is  computationally  very  expensive  for  flutter  calculations,  (2)  the  procedures  are  not 
easily  user-adaptive,  (3)  Structure  engineers  prefer  to  work  with  solutions  in  frequency-domain 
not  time-domain.  Industry  standard  methods  such  as  the  classical  AIC  method  is  lacking  in  the 
transonic  regime  where  a  reliable  high-level  CFD  method  such  as  CFL3D  is  badly  needed.  The 
proposed  frequency-domain  POD/ROM  method  originated  by  Hall/Dowell  has  a  great  potential 
to  be  such  transonic  AlC-like  method  which  has  the  essence  of  efficiency,  accuracy,  modularity 
and  could  assess  the  key  physics  for  nonlinear  flutter/LCO.  Its  future  applicability  is  far 
reaching  which  extends  to  areas  in  aeroelasticity,  aeroacoustics,  turbomachinery,  ASE,  MDO, 
etc.  On  the  other  hand,  the  ZAERO  aeroelastic  system  exclusively  adopts  the  AIC  methodology 
as  basis  for  the  unified  unsteady  aerodynamics.  Although  efficient,  ZAERO  is  incapable  to 
provide  nonlinear,  viscous  unsteady  aerodynamics  with  large  amplitude  oscillations.  Once 
developed,  the  transonic  frequency-domain  POD/ROM/HB  (Harmonic  Balance)  methodology 
along  with  its  Navier-Stokes  option,  would  be  an  integral  part  of  ZONA’s  ZAERO  aeroelastic 
system.  ZONA  envisions  that  the  production  code  of  the  frequency-domain  POD/ROM/HB 
method  should  become  the  future  industry  standard  for  routine  transonic  flutter/LCO  analysis. 

7.2  ZONA  Funded  Phase  III  for  Commercialization 

Stated  in  Section  7.0,  ZONA  will  fund  Duke  University  for  further  R&D  necessary  to  attain  a 
production  of  aeroelastic  software  for  commercialization.  To  do  so,  ZONA  will  utilize  its 
current  annual  license  income  fund  from  MSC/ZONA51  and  ZAERO  to  support  the  Phase  III 
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R&D  activity.  Industry  P-sites  will  be  set  up  with  Boeing/St.  Louis  and  Lockheed  Martin/Forth 
Worth  to  upgrade  the  code  as  a  viable  industry  production  software. 

7.3  Commercialization  of  the  Software  with  ZAERO 

When  the  frequency-domain  POD/ROM/HB  code  is  successfully  developed,  ZONA  will 
package  the  POD/ROM/HB/Navier-Stokes  methodology  (with  either  software  option)  into  a 
well-defined  commercial  software  product.  To  facilitate  its  sales,  ZONA  plans  to  couple  the 
code  with  the  unified  aerodynamic  module  within  ZAERO  and  jointly  market  for  both  codes. 

ZONA  has  currently  teamed  up  with  a  major  FEM  software  house  (MSC  Software)  for  ZONA51 
(in  MSC/NASTRAN)  and  ZAERO  licensing.  ZONA  currently  supports  over  120  users 
worldwide  for  the  MSC/NASTRAN  Aero  Option  II.  Released  in  April  1999,  ZAERO  has 
already  had  a  dozen  users  including  Lockheed  Martin/Vought  Systems,  Boeing/Commercial, 
DLR/Gottingen,  SDRC,  NASA/MSFC,  NASA/Dryden,  Scaled  Composites  Coleman  Aerospace, 
etc.  With  ZONA’s  extensive  aerospace  contractual/software  licensing  experience,  this  software 
product  is  expected  to  break  into  the  aerospace  market  quite  readily.  Prospective  customers 
include  ZONA’s  current  clientele,  other  aerospace  companies  and  DoD  organization. 

Other  commercialization  steps  include  advertising  through  magazines,  the  internet,  and  attending 
AIAA  conventions.  All  these  are  currently  being  pursued  by  ZONA  for  the  licensing  sales  of 
ZAERO. 
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CHAPTER  8 


CONCLUSIONS 


Phase  I  of  our  STTR  effort  has  concluded  successfully  with  a  number  of  achievements  attained 
and  a  clear  plan  developed  for  Phase  II  and  III  to  follow  upon  these  accomplishments.  The 
challenge  at  the  start  of  Phase  I  was  to  assess  the  viability  of  current  time  marching  CFD  codes 
for  transonic  flutter  and  limit  cycle  oscillation  (LCO)  analysis  and  design  and  also  to  develop 
computationally  more  efficient  and  physically  accurate  descriptions  of  aeroelastic  systems 
including  the  representation  of  the  unsteady  aerodynamic  forces.  This  has  been  done. 

Consistent  with  the  experience  and  observations  of  other  researchers  and  practioners,  we  have 
found  that  time  marching  CFD  codes  of  the  conventional  sort  are  simply  too  computationally 
intensive  to  be  viable  for  aeroelastic  modeling  in  an  industrial  application.  (See  the  discussion  in 
Chapter  5.)  Moreover  if  such  codes  were  to  be  used  for  the  development  of  actively  controlled 
smart  structures,  their  cost  would  be  prohibitive.  Our  studies  were  conducted  with  a  state-of-the- 
art  CFD  code,  CFL3D,  made  available  to  us  by  the  Aeroelasticity  Branch  of  the  NASA  Langley 
Research  Center.  We  are  most  appreciative  of  their  help  in  making  this  evaluation. 

Much  of  our  effort  in  Phase  I  has  therefore  been  devoted  to  the  development  of  alternative 
methods  using  the  concept  of  aerodynamic  modes,  either  eigenmodes  of  the  flow  field  or  the 
modes  one  can  construct  using  the  method  of  proper  orthogonal  decompostion  (POD).  Using  a 
few  of  the  most  dominant  modes,  Reduced  Order  Models  (ROM)  can  be  obtained  that  are  orders 
of  magnitude  smaller  in  terms  of  the  number  of  degrees  of  freedom  and  computational  time 
compare  to  the  original  CFD  model.  Such  models  are  simply  called  POD/ROM.  It  is 
emphasized  that  these  models  retain  the  essential  physical  modeling  from  the  CFD  code  from 
which  they  are  obtained,  e.g.  Euler  or  Navier-Stokes  equations. 

A  companion  development  has  been  the  study  of  large  shock  motions  using  the  method  of 
harmonic  balance  (HB)to  represent  the  dynamically  nonlinear  shock  wave  motion  and 
aerodynamic  forces  due  to  airfoil  (or  wing)  motion.  As  part  of  the  planned  work  for  Phase  II/III 
we  plan  to  combine  the  best  features  of  the  POD/ROM  and  HB  methods. 

Using  POD/ROM  and  HB  several  studies  were  undertaken  in  Phase  I  to  illustrate  the 
computational  efficiency  and  physical  modeling  capability  of  this  new  approach.  These  include 

•  Transonic  Limit  Cycle  Oscillation  Analysis  of  an  Airfoil  with  Control  Surface  Freeplay 
(Chapter  2) 

Here  the  POD/ROM  aerodynamic  forces  are  combined  with  a  nonlinear  structural  model  to 
investigate  LCO  due  to  freeplay. 

•  Nonlinear  Inviscid  Aerodynamic  Effects  on  Transonic  Divergence,  Flutter  and  Limit  Cycle 
Oscillations  (Chapter  3) 
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Here  the  harmonic  balance  (HB)  approach  is  used  to  model  the  aerodynamic  flow  with  large 
shock  wave  motions.  An  LCO  study  for  an  airfoil  oscillating  in  pitch  where  the  nonlinearity 
is  entirely  due  to  large  shock  motion  has  been  completed. 

•  Three  Dimensional  Transonic  Aeroelasticity  Using  Proper  Orthogonal  Decomposition  Based 
Reduced  Order  Models  (Chapter  4) 

For  small  shock  motions,  the  POD/ROM  methodology  has  been  extended  to  three 
dimensional  flows  over  wings  and  a  transonic  flutter  analysis  has  been  completed. 

Future  plans  are  described  in  more  detail  in  Chapter  6  and  will  include  the  addition  of  viscosity 
in  both  the  two-dimensional  and  three-dimensional  flow  POD/ROM  with  the  HB  method  to  be 
further  developed  to  account  for  large  (as  well  as  small)  shock  motions  and  also  for  viscous 
nonlinearities.  A  systematic  plan  and  roadmap  has  been  constructed  for  this  work  and  a 
commercializtion  strategy  developed  as  well  (Chapter  7). 
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INTRODUCTION 

The  principal  focus  of  this  paper  is  on  the  transonic  aeroe- 
lastic  behavior  of  an  airfoil  with  control  surface  freeplay  in¬ 
cluding  flutter  and  limit  cycle  oscillations.  In  most  of  this 
work,  we  will  assume  the  shock  motion  is  sufficiently  small 
such  that  it  is  (linearly)  proportional  to  the  airfoil  motion,  e.g. 
airfoil  motions  are  less  than  the  equivalent  of  one  degree  in 
angle  of  attack. 

Using  an  Euler/CFD-based  reduced  order  aerodynamic 
model,  a  thorough  study  of  the  flutter  boundary  with  Mach 
number  (M)  is  first  presented  in  the  absence  of  freeplay.  Par¬ 
ticularly  noteworthy  are  the  rapid  changes  of  flutter  modal 
content  in  the  transonic  range.  This  is  attributed  in  part  to 
the  rapid  changes  of  center  of  pressure  location  as  the  mean 
shock  position  changes  with  Mach  number.  These  changes  in 
the  modal  response  content  are  also  found  in  the  limit  cycle  os¬ 
cillations  (LCO)  which  are  encountered  when  control  surface 
freeplay  is  present.  Indeed  for  LCO,  the  modal  content  may 
change  at  a  fixed  Mach  number  when  the  dynamic  pressure  or 
flow  density  is  varied. 

Below  MO.80,  the  LCO  and  flutter  oscillations  are  qual¬ 
itatively  similar  to  those  found  at  low  Mach  number  where 
earlier  analyses  and  experiments  have  been  carried  out.  How¬ 
ever,  the  response  behavior  in  the  transonic  flow  regime  is  no¬ 
tably  different.  Of  special  interest  is  the  occurrence  of  flutter 
in  a  narrow  range  of  Mach  number  for  pitch  and  flap  (control 
surface)  dominated  motions.  Moreover,  beyond  a  certain  high 
transonic  Mach  number  (after  the  mean  shock  position  reaches 
the  trailing  edge  of  the  airfoil),  neither  flutter  nor  limit  cycle 
oscillation  occurs. 

SIGNIFICANCE  OF  LCO 

LCO  is  known  to  occur  on  various  operational  aerospace 
flight  vehicles.  This  has  been  a  source  of  serious  concern  since 
there  arc  no  analysis  techniques  available  that  have  predicted 
LCO  in  an  operational  aircraft.  There  have  been  some  semi- 
empirical  techniques  developed  to  correlate  with  LCO  that 
have  been  observed  in  flight,  and  these  are  useful  for  under¬ 
standing  the  LCO  that  has  occurred.  See  Reference  [1].  How- 
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ever  these  techniques  are  not  as  satisfactory  for  the  design  of  a 
new  vehicle  or  the  substantial  modification  of  an  existing  one, 
e.g.  new  stores  to  be  carried  by  an  aircraft. 

In  wind  tunnel  tests  of  flight  vehicle  prototypes,  LCO  has 
been  notably  absent  for  the  most  part.  This  is  perhaps  not  al¬ 
together  surprising  since  wind  tunnel  scale  models  have  been 
designed  based  upon  linear  aeroelastic  concepts.  Such  wind 
tunnel  models  and  tests  have  been  used  successfully  for  many 
years  to  predict  flutter  (the  onset  of  a  dynamic  linear  instabil¬ 
ity).  However  they  are  inherently  unable  to  predict  a  nonlinear 
phenomenon  such  as  LCO. 

In  this  regard,  it  should  be  noted  that  LCO  may  be  benefi¬ 
cial  as  well  as  detrimental.  Without  the  nonlinearities  that  lead 
to  LCO,  the  onset  of  flutter  may  lead  to  catastrophic  failure  of 
the  structure.  Hence  if  we  can  understand  and  predict  LCO, 
perhaps  we  can  take  advantage  of  these  nonlinearities  to  shape 
more  favorable  responses  of  the  aircraft  leading  to  enhanced 
safety  and  performance. 

TECHNICAL  DISCUSSION 
Sources  of  Nonlinearities 

The  principal  sources  of  the  nonlinearities  essential  to  the 
LCO  are  a  subject  of  current  debate  among  the  experts  in  the 
field.  The  candidate  sources  are  several: 

Fluid 

•  Shock  motions 

•  Separated  flow  motions 

Structure 

•  Free-play 

•  Geometric,  e.g.  a  nonlinear  relationship  between 
strain  and  displacement 

•  Material,  e.g  dry  friction  damping 

Also  there  is  a  further  distinction  between  a  static  versus  a  dy¬ 
namic  nonlinearity.  An  important  example  of  this  is  the  role 
of  a  shock  wave  in  the  fluid.  If  a  shock  is  present,  then  its 
creation  is  the  result  of  a  dynamic  nonlinear  process.  However 
once  a  steady  flow  is  established,  and  if  the  airfoil  motion  is 
sufficiently  small,  then  the  shock  motion  will  also  be  small  and 
proportional  to  the  airfoil  motion.  Hence  in  this  situation,  the 
shock  itself  represents  a  nonlinear  static  (time  independent) 
equilibrium  and  the  motion  may  be  treated  as  a  dynamically 
linear  perturbation  about  the  mean  shock  position.  In  most 
of  the  following  discussion,  we  assume  a  dynamically  linear 
model  of  the  shock  motion,  but  also  include  a  structural  (dy¬ 
namic)  nonlinearity,  i.e.  freeplay  in  the  connection  of  the  con¬ 
trol  surface  to  the  airfoil. 
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Figure  1:  Airfoil  with  Control  Surface  Configuration 


Figure  2:  Restoring  Moment  As  a  Function  of  Control  Surface 
or  Flap  Rotation  Angle 


Airfoil  with  Control  Surface  Freeplay  Model 

A  sketch  of  the  configuration  is  shown  in  Figure  (1).  It 
is  a  conventional  typical  section  model  except  that  the  spring 
that  attaches  the  control  surface  to  the  airfoil  has  a  nonlinear 
freeplay.  The  elastic  restoring  torque  or  moment  provided  by 
this  spring  is  shown  in  Figure  (2)  as  a  function  of  control  sur¬ 
face  or  flap  rotation  angle,  (3.  The  freeplay  angle  is  6 .  Note 
that  when  (3  is  less  than  S ,  there  is  zero  restoring  torque,  while 
for  (3  greater  than  6,  the  spring  stiffness  is  the  nominal  value 
in  the  absence  of  freeplay.  The  freeplay  may  be  thought  of 
as  creating  a  stiffness  or  (uncoupled)  natural  frequency  of  the 
spring  that  varies  as  a  function  of  flap  amplitude.  This  in¬ 
terpretation  is  shown  in  Figure  (3).  Here  the  flap  uncoupled 
natural  frequency  normalized  by  the  nominal  value  in  the  ab¬ 
sence  of  freeplay  is  shown  as  a  function  of  flap  amplitude,  /3> 
normalized  by  the  freeplay  angle,  6.  Note  that  given  a  certain 
flap  amplitude,  there  is  a  corresponding  ’’equivalent'’  flap  fre¬ 
quency.  Of  course  for  a  linear  system  the  control  surface  or 
flap  frequency  would  have  a  fixed  value  independent  of  flap 
amplitude. 


Figure  3:  Equivalent  Nonlinear  Flap  Frequency  Versus  Flap 
Rotation 


Thus  conceptually  and  computationally  one  may  proceed 
as  follows.  First,  one  determines  the  neutrally  stable  motions 
of  the  system  in  the  absence  of  freeplay  for  various  flap  fre¬ 
quencies  from  zero  to  the  nominal  value.  Then  one  determines 
the  corresponding  neutrally  stable  nonlinear  limit  cycle  mo¬ 
tion,  namely  the  flap  amplitude,  from  Figure  (3).  More  details 
of  this  procedure  will  be  given  later  including  the  construc¬ 
tion  of  the  linear  instability  or  flutter  boundary  for  this  system, 
which  serves  as  a  basis  for  determining  the  limit  cycle  oscilla¬ 
tions. 

This  approach  has  been  used  successfully  at  low  Mach 
number  and  the  theoretical  results  correlated  with  experi¬ 
ments.  See  References  [2]  and  [31.  To  extend  these  earlier 
calculations  to  transonic  flow  requires  a  substantially  more 
sophisticated,  but  still  computationally  efficient  aerodynamic 
model.  Fortunately  one  is  available  as  is  described  next. 

Computational  Fluid  Dynamic  (CFD)  Modeling  and 
its  Modal  Decomposition 

A  typical  CFD  model  is  very  large  in  terms  of  the  num¬ 
ber  of  equations  required  to  be  solved.  And  this  makes  such 
models  problematical  for  aeroelastic  (and  some  other)  analy¬ 
ses.  For  example,  the  CFD  model  used  in  the  present  work  is 
based  upon  the  Euler  equations  of  fluid  mechanics  and  has  a 
spatial  grid  of  65  X  97  (6035)  mesh  points.  At  each  grid  point 
there  are  four  fluid  variables  to  be  determined.  Thus  the  CFD 
model  per  se  has  about  25,000  unknowns  to  be  determined  by 
solving  25,000  equations.  This  is  a  doable  task  if  the  struc¬ 
tural  motion  is  known.  However  if  this  CFD  model  is  to  be 
combined  with  a  set  of  structural  equations  of  motion,  and  so¬ 
lutions  are  to  be  found  for  many  combination  of  structural  and 
fluid  parameters,  then  the  calculation  using  the  original  CFD 
model  quickly  gets  out  of  hand.  Thus  the  search  for  an  alter¬ 
native  approach. 
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(a)  Baseline  (b)  Grid  Motion 


Figure  4:  NACA  0012  Grid:  65  X  97  Computational  Nodes, 
Outer  Boundary  Radlus=15c 


Several  semi-empirical  methods  have  been  derived  to  ad¬ 
dress  the  computational  feasibility  issue.  These  as  well  as  the 
method  to  be  described  and  used  here  are  discussed  in  more 
depth  in  Reference  [4],  Among  these  methods  are  variations 
on  the  notion  of  a  Pade  approximant. 

The  method  used  here  is  based  upon  the  observation  that 
virtually  all  CFD  models  can  be  thought  of  as  having  a  modal 
composition.  The  simplest  conceptual  set  of  modes  is  perhaps 
the  fluid  or  aerodynamic  eigenmodes  of  the  CFD  model,  and 
these  modes  have  been  used  successfully  in  creating  reduced 
order  models  (ROM)  that  are  computationally  and  conceptu¬ 
ally  attractive.  See  the  discussion  in  Reference  [4]  and  the 
forthcoming  Reference  [5]. 

However  it  turns  out  that  determining  the  aerodynamic 
eigenmodes  of  a  large  CFD  model  is  itself  a  challenging 
task.  Hence  a  method  called  Proper  Orthogonal  Decompo¬ 
sition  (POD)  is  employed  here  and  in  Reference  [4].  The  first 
use  of  this  method  in  an  Euler  based  aerodynamic  context  was 
by  Romanowski.  Again  see  the  discussion  of  the  literature  in 
Reference  [4].  In  this  method,  by  using  a  time  history  or  fre¬ 
quency  response  of  the  CFD  model  to  a  known  structural  in¬ 
put,  one  may  construct  a  small  set  of  modes  or  basis  functions, 
typically  less  than  100.  Using  these  modes,  one  can  reconsti¬ 
tute  and  very  substantially  reduce  the  size  of  the  CFD  model 
with  essentially  no  loss  in  accuracy  or  physical  content.  This 
is  the  approach  used  here. 

For  the  present  analysis,  63  POD  modes  are  found  from  the 
frequency  responses  (aerodynamic  transfer  functions)  in  flap, 
pitch  and  plunge  respectively  at  21  frequencies  at  each  Mach 
studied  using  the  original  CFD  model.  Based  upon  previous 
experience,  one  might  use  an  even  smaller  number  of  aero¬ 
dynamic  modes  than  this.  However,  even  with  this  generous 
number  of  modes,  the  computations  described  below  were  all 
done  in  a  few  days. 

The  computational  grid  used  for  the  CFD  model  is  shown 
in  Figure  (4),  and  Figure  (5)  shows  the  chordwise  steady  flow 
pressure  distributions  for  several  Mach  numbers.  Note  the 
presence  of  the  shock  at  M=0.80,  0.85,  and  0.90. 

Linear  Instability  (Flutter) 

First,  consider  the  flutter  behavior  for  this  system  in  the 
absence  of  freeplay.  The  stability  of  this  system  was  assessed 


Figure  5:  Computed  NACA  0012  Steady  Flow  Surface  Pres¬ 
sure  Distributions:  a o  =  0.0  (deg) 


Structural  Parameters 

Half  Chord 

b  = 

0.127m 

Sectional  Mass* 

m  = 

1.567  kg/m 

Air  Density 

Poc  — 

1.225  kg/m' 

Pitch  Axis  Location  (e/b) 

a  = 

—0.5 

Flap  Hinge  Location  (ef/b) 

af  = 

0.5 

Mass  Ratio,  (m/p^nb) 

P  = 

25.24 

Static  Unbalance,  ( Sa/mb ) 

2q 

0.4316 

Radius  of  Gyration,  ( Ia/mb 2) 

r2  = 

0.5331 

Flap  Static  Unbalance,  (Sp  /mb) 

xp  = 

0.01985 

Flap  Radius  of  Gyration,  (Ip  / mb2) 

Tl  = 

0.01292 

Uncoupled  Frequencies 

uh  =  4.59  Hz 

ua  =  6.04  Hz 

up  =  12.54  Hz 

Coupled  Frequencies 

a;*  =  4.45  Hz 

wQ  =  9.21  Hz 

up  =  19.44  Hz 

*Plunge  inertia  corrected  for  experimental  support  mass 
M/m  =  2.166 

Table  1 :  Structural  Parameters  for  NACA  0012  Airfoil  with  Con¬ 
trol  Surface 

by  constructing  a  root  locus  (migration  of  the  true  aeroelastic 
eigenvalues)  as  a  function  of  the  nondimensional  airspeed  or 
dynamic  pressure  for  each  Mach  number.  The  usual  structural 
and  flow  parameters  are  defined  in  Table  1 .  These  are  typi¬ 
cal  and  correspond  to  the  theoretical  and  experimental  model 
studies  in  References  [1]  and  [2]. 

A  representative  root  locus  result  is  shown  in  Figure  (6)  for 
M=0.80.  Root  locus  results  are  available  for  all  Mach  num- 
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Real  Frequency  Ratio,  Re(X/coa) 


Figure  6:  Aeroelastic  System  Eigenvalue  Root-Loci:  M  =  0.80, 
a0  =  0.0  (deg) 

bers  used  to  construct  the  flutter  boundary  which  is  shown  in 
Figure  (7).  Consider  first  Figure  (6).  The  eigenvalues  in  this 
root  locus  are  those  of  the  coupled  fluid/structural  (aeroelastic) 
system.  However  the  roots  that  originate  as  structural  modes 
at  low  values  of  velocity  or  dynamic  pressure  are  readily  iden¬ 
tified  and  labeled  in  the  figure.  The  vertical  axis  is  the  imag¬ 
inary  part  of  the  eigenvalue  or  frequency,  and  the  horizontal 
axis  displays  the  real  part  or  rate  of  growth  (if  positive)  or  de¬ 
cay  (if  negative)  of  the  oscillations  associated  with  each  root 
or  eigenvalue. 

Note  especially  that  the  plunge”  aeroelastic  mode  has  a 
root  that  for  low  ’’gain”  (or  flow  velocity  or  dynamic  pres¬ 
sure)  moves  to  the  left  and  becomes  more  stable.  But  then 
as  the  flow  velocity  increases,  it  reverses  direction  and  moves 
into  me  right  half  plane  becoming  unstable.  And  then  at  even 
higher  velocities,  it  moves  back  into  the  left  hand  plane  and 
becomes  stable  again.  However  by  then,  the  pitch  mode  has 
moved  into  the  right  half  plane  and  become  unstable.  Hence 
the  aeroelastic  system  remains  unstable  once  the  plunge  mode 
becomes  again  stable  at  this  velocity  and  Mach  number. 

All  the  other  roots  in  this  figure  which  appear  to  originate 
from  near  the  origin  are  essentially  aerodynamic  roots,  and 
these  roots  all  move  off  into  the  left  half  plane  indicating  they 
are  always  stable  and  increasingly  so  as  the  dynamic  pressure 
increases.  On  occasion  an  aerodynamic  root  may  become  un¬ 
stable  however,  though  not  for  the  parameters  studied  in  this 
work. 

As  the  Mach  number  becomes  higher,  the  most  critical  root 
may  change.  For  example,  at  M=0.85  the  pitch  root  becomes 
unstable  first,  and  for  M=0.90,  it  is  the  flap  root.  At  yet  higher 
Mach  numbers,  no  roots  become  unstable.  For  brevity  these 
other  root  loci  are  omitted  here. 

Taking  all  of  this  information  from  the  root  loci  at  various 
Mach  numbers,  the  flutter  boundary  trend  with  Mach  num- 


Mach  Number,  M 


(a)  Flutter  Speed 


Mach  Number,  M 

(b)  Flutter  Frequency 


Figure  7:  Mach  Number  Flutter  Trend:  ao  =  0.0  (deg) 


ber  can  be  determined,  and  is  presented  in  Figure  (7a).  There 
are  several  interesting  features  to  this  flutter  boundary.  Up  to 
M=0.80,  the  root-loci  are  rather  similar,  and  it  is  always  the 
plunge  root  that  is  critical  for  flutter.  Starting  at  M=0.80,  the 
pitch  root  also  shows  instability,  and  at  M=0.825  and  0.85,  it 
is  most  critical  for  flutter.  At  M=0.875  and  0.90,  the  flap  mode 
is  most  critical  for  flutter,  and  for  M=0.925  to  at  least  M=1 .1, 
no  flutter  is  observed  for  a  non-dimensional  flow  velocity  up 
to  at  least  one.  The  corresponding  frequencies  of  the  flutter 
oscillations  are  shown  in  Figure  (7b).  Note  that  in  Figure  (7) 
when  two  data  points  are  shown  for  say  the  plunge  root  at  a 
fixed  Mach  number,  the  lower  velocity  point  is  when  flutter 
begins,  and  the  higher  velocity  point  is  when  the  root  returns 
to  the  stable  left  half  plane  and  flutter  ceases  in  that  root. 

Note  also  the  narrow  range  of  Mach  number  where  the 
change  in  flutter  mode  occurs.  Results  of  this  type  have  been 
observed  in  experiments  where  they  are  called  ’’chimneys”. 
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See  Reference  [6]. 

Limit  Cycle  Oscillations 

Now  the  ffeeplay  is  added  to  the  model  and  thus  LCO  may 
occur  As  is  perhaps  obvious  from  physical  intuition,  when 
ffeeplay  is  added,  the  stiffness  of  the  control  surface  ffeeplay 
is  reduced  for  small  motions!  Hence  one  expects  limit  cycle 
oscillations  to  occur  below  the  flutter  boundary,  i.e.  at  flow 
velocities  less  than  those  shown  in  Figure  (7a).  Indeed  a  few 
moments  of  reflection  may  lead  one  to  expect  that  once  the 
linear  flutter  boundary  shown  in  Figure  (7a)  is  exceeded,  then 
exponentially  explosive  flutter  will  occur  when  the  nonlinear¬ 
ity  is  due  to  freeplay.  That  is  found  to  be  the  case  as  shown  by 
the  present  analysis  and  also  by  the  analysis  and  experiments 
of  References  [2]  and  [3]. 

To  compute  the  LCO  one  proceeds  as  follows.  For  all  other 
parameters  fixed  including  Mach  number,  the  flutter  veloc¬ 
ity  is  determined  from  a  linear  flutter  analysis  as  a  function 
of  (uncoupled)  flap  natural  ffequency  (or  equivalently  spring 
stiffness).  Since  the  flap  ffequency  (or  stiffness)  is  known 
as  a  function  of  flap  amplitude,  see  Figure  (3  (or  (2)),  we 
immediately  can  determine  from  these  two  results  (by  cross- 
plotting),  the  flap  amplitude  for  neutrally  stable  nonlinear  mo¬ 
tions  (LCO)  as  a  function  of  flow  velocity.  As  an  example 
for  M=0.80,  see  Figure  (8).  In  Figure  (8a),  the  ’’flutter”  or 
flow  velocity  at  which  neutrally  stable  oscillations  may  occur 
is  shown  as  a  function  of  flap  ffequency.  Using  Figure  (3), 
which  shows  the  nonlinear  dependence  of  flap  ffequency  as 
a  function  of  flap  amplitude,  Figure  (8b)  may  be  constructed. 
This  shows  LCO  amplitude  as  a  function  of  flow  velocity.  The 
corresponding  ’’flutter”  ffequency  and  the  LCO  frequency  are 
shown  in  Figures  (8c)  and  (8d)  respectively 

A  comment  on  this  method  is  appropriate  here.  The  reader 
familiar  with  the  harmonic  balance  approach  will  recognize 
that  a  single  harmonic  approximation  has  been  used  to  de¬ 
scribe  the  freeplay  nonlinearity.  This  approach  is  more  fully 
described  in  the  present  context  in  Reference  [3].  Of  course 
this  is  a  classical  procedure  for  nonlinear  systems;  though  it  is 
not  often  used  for  systems  with  as  many  degrees  of  freedom  as 
the  present  one.  The  results  of  Refei  nee  [3]  confirm  that  the 
harmonic  balance  approach  gives  results  in  good  agreement 
with  those  obtained  from  time  marching  solutions  that  include 
all  harmonics,  as  well  as  the  results  from  experiments. 

The  results  of  Figure  (8)  have  several  interesting  features. 
First  of  all,  the  limit  cycle  amplitude  is  normalized  by  the 
freeplay  angle,  5.  The  theoiy  predicts  and  experiments  agree, 
see  Reference  [3],  that  when  the  results  are  normalized  in  this 
manner,  they  are  universal.  That  is,  the  limit  cycle  amplitude 
is  proportional  to  the  freeplay  angle.  Secondly,  the  lowest  ve¬ 
locity  at  which  LCO  may  occur  corresponds  to  the  minimum 
flutter  velocity  that  occurs  at  a  certain  flap  frequency  See 
Figure  (8a)  and  then  compare  the  lowest  velocity  for  LCO  in 
Figure  (8b).  Strictly  speaking,  a  finite  disturbance  is  required 
to  generate  LCO  at  this  lowest  velocity  and  for  a  small  ve¬ 
locity  range  thereafter.  LCO’s  for  any  disturbance,  no  matter 
how  small,  will  only  occur  when  the  flutter  velocity  for  a  flap 
ffequency  of  zero  is  exceeded.  Again  compare  Figures  (8a) 
and  (8b).  The  unstable  LCO’s,  which  are  shown  along  with 
the  stable  LCO’s  (those  that  are  observed  in  an  experiment), 
provide  a  measure  of  the  level  of  disturbance  required  to  ini¬ 


tiate  the  LCO  at  the  lower  flow  velocities.  In  practice  such 
disturbances  are  usually  present  in  wind  tunnel  experiments 
and  flight  operations. 

The  results  of  Figure  (8)  are  typical  until  one  reaches  the 
higher  transonic  Mach  numbers  where  linear  theory  predicts 
flutter  will  cease.  At  the  highest  Mach  number  considered  here 
where  flutter  and  LCO  may  occur,  MK).90,  the  LCO  has  a 
somewhat  different  character.  See  Figure  (9).  Again  the  LCO 
is  first  encountered  at  the  minimum  velocity  at  which  flutter 
will  occur  over  the  range  of  flap  frequencies.  But  now  the  cor¬ 
responding  flap  frequency  is  zero.  See  Figure  (9a).  Moreover, 
when  the  flow  velocity  increases  to  higher  values,  there  are 
two  stable  limit  cycles.  The  nature  of  the  disturbances  to  the 
system  would  determine  which  of  these  two  LCO  would  be 
obsen/ed  in  a  wind  tunnel  experiment  or  in  flight.  Note  that 
the  LCO  branch  with  the  larger  amplitude  will  again  move  to 
amplitudes  with  very  large  values  (to  infinity  according  to  the 
present  theoretical  model)  when  the  flow  velocity  approaches 
the  linear  flutter  velocity. 

CONCLUSIONS 

The  transonic  flutter  and  limit  cycle  oscillations  of  an  air¬ 
foil  with  control  surface  freeplay  have  been  determined  using 
a  new  aerodynamic  modeling  technique  that  provides  greater 
physical  insight  and  understanding  by  tracing  the  true  root  lo¬ 
cus  of  the  corresponding  linear  aeroelastic  system.  This  in 
turn  enables  a  very  computationally  efficient  harmonic  bal¬ 
ance  technique  to  be  used  in  determining  the  nonlinear  limit 
cycle  oscillations. 

New  physical  insights  gained  include  the  rapid  change  in 
flutter  mode  that  occurs  in  the  transonic  Mach  number  range. 
This  phenomenon  has  been  observed  in  experiments,  but  has 
not  been  previously  predicted  theoretically.  With  respect  to 
LCO,  these  are  the  first  results  available  in  the  transonic  range 
for  the  configuration  studied.  The  model  also  predict?  signifi¬ 
cant  changes  in  LCO  behavior  as  a  function  of  Mach  number. 
However,  these  are  as  yet  unconfirmed  by  experiments.  Up  to 
high  subsonic  Mach  numbers,  the  flutter  and  LCO  results  are 
similar  to  those  previously  found  at  low  Mach  numbers. 

CURRENT  WORK 

Current  work  is  directed  toward  addressing  LCO  arising 
from  large  amplitude  shock  motions  (and  even  in  the  absence 
of  any  structural  nonlinearity  such  as  freeplay).  Preliminary 
results  have  now  been  obtained,  and  we  have  determined  LCO 
triggered  by  a  single  degree  of  freedom  flutter  in  pitch  in  the 
transonic  range.  For  the  particular  configuration  studied,  flut¬ 
ter  may  occur  below  the  linear  flutter  velocity  if  the  distur¬ 
bance  is  sufficiently  large,  say  on  the  order  of  3  degrees  of 
pitch.  These  results  are  described  in  a  separate  abstract  that  is 
being  submitted  for  presentation  at  the  SDM  Conference. 
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(b)  LCO  Flap  Rotation 
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INTRODUCTION 

Limit  cycle  oscillations  (LCO)  in  aeroelastic  systems  ap¬ 
pear  to  be  more  prevalent  in  transonic  flow  than  in  subsonic 
flow.  Hence  it  has  been  thought  that  at  least  for  some  configu¬ 
rations  the  source  of  the  nonlinearity  that  leads  to  LCO  is  in  the 
aerodynamic  flow.  Of  course  nonlinear  structural  mechanisms 
can  also  lead  to  LCO  whether  the  flow  is  transonic  or  not.  And 
there  have  been  wind  tunnel  experiments  where  the  test  model 
was  designed  to  exhibit  LCO  due  to  a  structural  nonlinearity, 
and  such  test  results  have  been  successfully  correlated  with 
analysis.  See  References  [1]  and  [2].  However,  the  present 
understanding  of  LCO  induced  by  aerodynamic  nonlinearities 
is  less  complete  and  as  yet  no  systematic  quantitative  correla¬ 
tion  between  theory  and  experiment  has  been  achieved. 

This  is  perhaps  a  meaningful  measure  of  the  greater  dif¬ 
ficulty  in  modeling  aerodynamic  non  linearities,  both  theoreti¬ 
cally  and  experimentally,  compared  to  modeling  nonlinearities 
in  a  structure. 

One  of  the  advantages  of  studying  theoretical  models  is 
that  each  of  the  several  possible  physical  phenomena  that  may 
lead  to  LCO  can  be  studied  separately.  In  this  paper,  we  con¬ 
sider  the  effects  of  nonlinearities  arising  from  inviscid  tran¬ 
sonic  aerodynamics.  The  principal  physical  effect  of  interest 
is  the  relatively  large  motion  of  the  shock  wave  as  the  am¬ 
plitude  of  say  the  pitch  motion  of  the  airfoil  becomes  suffi¬ 
ciently  large.  This  in  turn  leads  to  a  movement  of  the  center 
of  pressure  with  amplitude.  Hence  one  expects  to  see  an  effect 
of  amplitude  on  the  neutrally  stable  motions  that  may  occur. 
Moreover  this  may  lead  to  limit  cycle  motions  rather  than  the 
catastrophic  exponentially  growing  oscillations  predicted  by 
time  linearized  models.  The  latter  models  capture  the  effect 
of  the  mean  position  of  the  shock  and  small  shock  motions 
about  this  mean  position  by  assuming  the  shock  motion  is  dy¬ 
namically  linear,  i.e.  the  shock  motion  is  proportional  to  the 
airfoil  motion.  This  is  not  true  for  dynamically  nonlinear  aero¬ 
dynamic  models  that  allow  for  larger  and  more  general  shock 
motions  including  the  possible  appearance  and  disappearance 
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of  a  shock  during  a  cycle  of  airfoil  motion.  The  latter  is  our 
concern  here. 

TECHNICAL  DISCUSSION 

In  this  paper,  we  will  consider  two  distinct  aeroelastic  phe¬ 
nomena,  divergence  and  flutter,  and  their  associated  limit  cy¬ 
cle  oscillations.  To  keep  the  discussion  focussed  on  the  fun¬ 
damental  physical  phenomena,  and  to  ease  the  interpretation 
of  the  inherently  complex  phenomena,  only  a  single  structural 
degree  of  freedom  will  be  studied.  However  the  aerodynamic 
modei  will  be  a  state-of-the-art  computational  fluid  dynamics 
(CFD)  based  upon  the  Euler  equations  of  nonlinear,  rotational 
inviscid  aerodynamic  theory.  The  aerodynamic  model  and  its 
spatial  discretization  will  be  discussed  in  the  full  paper. 

Here  we  emphasize  that  the  solution  technique  is  for  a 
large  system  of  ordinary  differential  equations  in  time,  which 
represents  the  time  variation  of  the  fluid  unknowns  at  each  spa¬ 
tial  grid  point  in  the  CFD  model.  The  unknowns  are  four  in 
number  at  each  grid  point  for  a  two-dimensional  Euler  flow 
and,  for  example,  could  be  density,  the  two  scalar  components 
of  momentum,  and  the  total  energy  at  each  grid  point.  The 
present  CFD  model  has  about  17,000  total  flow  variable  un¬ 
knowns,  and  therefore  an  efficient  solution  method  is  impera¬ 
tive  to  carry  out  the  studies  reported  here. 

Harmonic  Balance  Solution  in  the  Frequency 
Domain 

The  pioneering  work  of  Ueda  and  Dowell  [3]  and  Lan  and 
his  colleagues  [4]  should  be  recalled.  Ueda  and  Dowell  used  a 
describing  function  technique  whereby  the  dominant  harmonic 
was  extracted  from  a  time  marching  CFD  model,  LTRAN2, 
using  both  indicial  and  harmonic  motions  of  the  airfoil.  They 
considered  a  two  degree  of  freedom  typical  airfoil  section.  Lan 
et.  al.  used  the  method  of  harmonic  balance  to  study  the  un¬ 
steady  transonic  aerodynamics  for  flutter  and  limit  cycle  oscil¬ 
lation  prediction.  In  their  work,  they  used  the  transonic  small 
disturbance  potential  flow  model,  as  did  Ueda  and  Dowell,  and 
only  considered  a  single  harmonic.  In  the  present  work,  we 
employ  the  Euler  equations  of  fluid  dynamics  and  also  retain 
multiple  harmonics  in  the  aerodynamic  model.  It  is  found  that 
using  several  harmonics  improves  the  theoretical  prediction  of 
the  aerodynamic  forces.  However  in  the  aeroelastic  analysis, 
when  the  fluid  and  structural  models  are  coupled,  only  a  single 
harmonic  is  used.  The  effects  of  higher  harmonics  on  this  sin¬ 
gle  harmonic  are  retained  as  they  are  found  to  be  significant  in 
the  fluid  model. 
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(a)  Close-up  (b)  Overall 


Figure  1:  NACA  64A010A  Computational  Grid 

The  Aeroelastic  System  and  Its  Solution 

The  structural  equation  of  motion  is  a  simple  single  degree 
of  freedom  model  in  pitch.  See  Figure  (1)  for  a  depiction  of 
the  airfoil  and  the  CFD  grid  used  in  the  numerical  calcula¬ 
tions.  By  carefully  selecting  the  pitch  axis  and  mass  ratio,  we 
can  insure  that  the  system  will  either  undergo  classical  linear 
aeroelastic  divergence  or  flutter.  Divergence  can  occur  when 
the  aerodynamic  "negative”  stiffness  overcomes  the  structural 
stiffness,  while  flutter  may  occur  when  the  aerodynamic  neg¬ 
ative  damping  overcomes  the  structural  damping.  As  will  be 
shown,  each  of  these  classical  linear  aeroelastic  phenomena 
has  a  distinctively  different  limit  cycle  or  nonlinear  behavior. 

The  Mach  number  for  these  studies  is  M=0.80  and  a  NACA 
64A010A  airfoil  is  considered.  The  NACA  64A010A  is 
a  symmetric  (10.6%  thickness  ratio)  variant  of  the  "Arnes” 
AGARD  156  benchmark  section.  The  elastic  axis  is  consid¬ 
ered  at  the  mid-chord.  Employed  for  the  CFD  calculations 
is  an  “0”-type  computational  mesh  with  65  x  65  radial  and 
circumferential  nodes  that  has  an  outer  boundary  radius  of  10 
chord  lengths.  The  computed  static  pressure  distribution  for  an 
angle  of  attack  of  0.0  and  5.0  degrees  is  shown  in  Figure  (2). 
Note  that  at  5.0  degrees,  the  upper  surface  shock  wave  has 
moved  rearward  and  increased  in  strength,  and  on  the  lower 
surface,  the  shock  has  essentially  disappeared.  The  center  of 
pressure  (zcp)  as  a  function  of  static  angle  of  attack  is  shown 
Figure  (3)  where  it  is  seen  the  center  of  pressure  moves  from 
32%  chord  to  40%  chord  as  the  angle  of  attack  varies  from  0.0 
to  5.0  degrees. 

Linear  and  Nonlinear  Divergence 

This  is  perhaps  the  simpler  of  the  two  phenomena  since 
by  definition  it  is  time  independent,  i.e.  we  are  dealing  with 
a  static  linear  instability  and  its  nonlinear  counterpart.  In  this 
case,  the  structural  equation  of  motion  becomes  an  equation  of 
static  equilibrium.  And  for  the  aerodynamic  model,  we  only 
need  to  determine  the  lift  and  moment  about  some  appropri¬ 
ate  axis  as  a  function  of  angle  of  attack.  For  small  angle  of 
attack,  we  will  recover  the  classical  linear  aeroelastic  diver¬ 
gence  phenomena.  But  the  question  is,  what  are  the  effects  of 
the  nonlinearity? 

The  equation  of  static  equilibrium  simply  equates  the  aero¬ 
dynamic  and  elastic  restoring  moments.  Namely, 

K &  —  QcoC  Cma  ( ^ ) 


Figure  2:  Steady  Flow  Surface  Pressure  Distributions:  NACA 
64A010A  Airfoil  Section,  M  =  0.80 


Fipure  3:  Center  of  Pressure  Variation  with  Angle  of  Attack: 
NACA  64A010A  Airfoil  Section,  M  =  0.80 

Defining  a  non-dimensional  dynamic  pressure,  A,  Equation  (1) 
may  be  rewritten  as 

A  =  a/cmQ  (2) 

where  A  is  given  by  A  =  qoo^/Kco  The  angle  of  attack  may 
have  an  initial  angle,  a0,  which  is  prescribed,  and  also  an  ad¬ 
ditional  angle  due  to  the  torsional  twist  of  the  elastic  spring, 
a. 

Now  for  a  linear  aeroelastic  model,  the  aerodynamic  mo¬ 
ment  coefficient  is  simply  proportional  to  the  angle  of  attack. 
Thus  for  no  initial  angle  of  attack,  the  classical  linear  diver¬ 
gence  dynamic  pressure  is  given  by  Equation  (2)  where  A  is 
now  a  fixed  number. 

To  extend  this  study  of  divergence  into  the  nonlinear  range, 
we  recognize  that  now  the  aerodynamic  coefficient  is  a  nonlin¬ 
ear  function  of  angle  of  attack.  For  zero  initial  angle  of  attack, 
we  may  determine  the  twist  of  the  torsional  spring,  and  its  de- 
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Nondimensional  Dynamic  Pressure,  X=q_c2/Ka  (radians) 

Figure  4:  Divergence  and  Post-Divergence  of  an  Airfoil  In¬ 
cluding  Transonic  Nonlinear  Inviscid  Aerodynamics:  NACA 
64A010A  Airfoil  Section,  M  =  0.80,  Elastic  Axis  Location, 
a  =  e/b  =  0.0 


pendence  on  A,  by  specifying  the  twist  angle  in  Equation  (2) 
and  then  solving  for  A.  This  procedure  is  readily  extended  to 
the  case  with  an  initial  angle  of  attack. 

Qualitatively  one  can  anticipate  the  effect  of  the  aerody¬ 
namic  nonlinearity  by  examining  the  aerodynamic  moment 
variation  with  angle  of  attack.  A  necessary  condition  for  di¬ 
vergence  to  occur  is  that  the  aerodynamic  moment  be  posi¬ 
tive  in  the  same  direction  as  the  twist  angle.  Moreover,  if  the 
nonlinear  aerodynamic  model  predicts  a  moment  less  in  mag¬ 
nitude  than  that  predicted  by  linear  aerodynamic  theory,  the 
effect  of  the  nonlinearity  will  be  to  stabilize  the  divergence. 
And  vice  versa  if  the  nonlinear  theory  predicts  an  increase  in 
aerodynamic  moment  over  that  given  by  linear  theory.  Hence 
by  examining  the  slope  of  the  moment  vs.  angle  of  attack 
curve  with  increasing  angle  of  attack,  we  will  know  whether 
the  effect  of  the  nonlinearity  is  favorable  or  unfavorable. 

In  the  example  below,  the  effect  is  favorable.  That  is,  once 
the  divergence  dynamic  pressure  for  a  small  angle  of  attack 
is  exceeded  (this  is  the  classical  linear  aeroelastic  divergence 
dynamic  pressure),  then  the  angle  of  twist  of  the  pitch  spring 
remains  finite  and  smoothly  increases  from  zero  beyond  the 
divergence  dynamic  pressure.  See  Figure  (4)  where  the  angle 
of  twist  is  plotted  vs  the  non-dimensional  dynamic  pressure. 
Also  shown  are  results  with  an  initial  angle  of  attack.  In  this 
latter  case,  there  is  some  twist  over  the  full  range  of  dynamic 
pressure.  Indeed  even  if  the  initial  angle  of  attack  is  only  a  few 
degrees,  it  would  be  difficult  to  detect  the  classical  divergence 
dynamic  pressure  experimentally  for  this  example.  For  readers 
who  have  studied  buckling  of  systems  in  the  presence  of  im¬ 
perfections  (e.g.  beams,  plates  or  shells  with  initial  curvature), 
this  behavior  will  be  familiar. 

In  this  example,  recall  the  center  of  pressure  moves  from 
32%  chord  at  low  angles  of  attack  to  40%  chord  at  5.0  degrees 
angle  of  attack.  This  is  the  principal  reason  for  the  stabiliz¬ 
ing  effect  of  nonlinear  aerodynamics  on  the  post-divergence 
condition. 

Had  the  change  of  the  slope  of  the  aerodynamic  moment 
curve  been  in  the  opposite  direction,  then  the  angle  of  twist 


vs.  dynamic  pressure  curve  would  have  bent  the  other  way. 
That  is,  for  dynamic  pressures  below  the  classical  divergence 
dynamic  pressure,  there  would  be  non-trivial  (non-zero)  twist 
angles  that  represent  possible  static  nonlinear  equilibrium  so¬ 
lutions.  Intuitively  one  recognizes  that  these  latter  solutions 
would  themselves  be  unstable,  i.e.  such  results  would  be  inter¬ 
preted  physically  as  the  magnitude  of  the  disturbance  required 
to  generate  non-trivial  twist  at  dynamic  pressures  below  the 
classical  divergence  dynamic  pressure.  In  our  studies  to  date, 
only  the  stable  nonlinear  effect  has  been  observed  for  stati¬ 
cally  divergent  systems.  However,  this  is  not  to  say  that  unsta¬ 
ble  nonlinear  divergence  systems  may  not  be  encountered  for 
some  other  parameter  combinations. 

Of  course,  divergence  is  a  very  special  case  of  nonlin¬ 
ear  aeroelasticity  as  it  is  for  linear  aeroelasticity,  because  the 
frequency  of  oscillation  is  zero  when  divergence  and  post¬ 
divergence  occurs.  Thus  we  now  turn  to  an  oscillatory  case. 

Flutter  and  Associated  LCO 

Now  consider  single-degree-of-freedom  flutter  in  pitch. 
Here  the  classical  flutter  arises  from  a  negative  damping  in 
the  aerodynamic  moment  beyond  a  certain  reduced  frequency. 
However  the  reduced  frequency  at  which  the  aerodynamic 
damping  moment  becomes  negative  increases  as  the  angle  of 
pitch  oscillation  increases.  Hence  the  reduced  velocity  de¬ 
creases  as  the  angle  of  pitch  increases,  which  suggests  that 
this  will  lead  to  an  unstable  LCO  as  indeed  it  does. 

In  the  example  considered,  we  have  moved  the  elastic  axis 
to  20%  chord  to  preclude  divergence  and  to  induce  flutter. 

It  should  be  emphasized  that  in  the  present  analysis,  we 
are  using  a  single  harmonic  to  represent  the  pitch  oscillation. 
However  in  the  calculation  of  the  aerodynamic  moment,  we 
include  up  to  three  harmonics  to  determine  the  effect  of  higher 
harmonics  on  the  first  harmonic  of  the  aerodynamic  moment. 
It  turns  out  that  the  effect  of  the  third  harmonic  is  negligible. 
Indeed,  if  one  only  retains  a  single  harmonic  in  the  aerody¬ 
namic  analysis,  the  results  are  qualitatively  correct  and  have 
fair  quantitative  accuracy. 

Results  for  the  first  harmonic  for  the  lift  and  moment  about 
the  pitch  or  elastic  axis  are  shown  in  Figure  (5).  These  re¬ 
sults  are  for  two  harmonics  retained  in  the  aerodynamic  anal¬ 
ysis.  Note  that  the  results  at  a  reduced  frequency  of  zero  were 
those  used  in  the  divergence  analysis  discussed  previously.  Of 
course,  a  transformation  of  the  pitch  axis  is  used  for  the  diver¬ 
gence  analysis. 

With  the  real  and  imaginary  parts  of  the  aerodynamic  mo¬ 
ment  taken  from  Figure  (5),  and  using  the  usual  pitch  equation 
of  motion, 

Ia  (q  +  2<at jaa  +  u>la)  =  q^Cmc  (3) 

where  =  Ka/Ia,  we  can  convert  this  equation  into  the 
frequency  domain,  and  nondimensional ize  and  separate  it  into 
real  and  imaginary  parts.  With  some  rearrangement,  these  two 
equations  can  be  written  as 
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(a)  Real  Unsteady  Lift  (b)  Imaginary  Unsteady  Lift 


Reduced  Frequency,  k=coc/U.  Reduced  Frequency,  k=coc/U„ 


(c)  Real  Unsteady  Moment 


(d)  Imaginary  Unsteady  Moment 


Figure  5:  Unsteady  Lift  and  Moment  for  Various  Pitch  Amplitudes:  NACA  64A010A  Airfoil  Section,  M  =  0.80,  qo  =  0.0  (deg), 
Elastic  Axis  Location,  a  =  e/b  =  -0.6,  Two  Harmonics  Employed  in  Harmonic  Balance  Expansion 
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where  a  bar  over  the  aerodynamic  coefficient  and  angle  of 
twist  denotes  the  amplitude  of  the  first  harmonic,  and  Re  and 
Im  denote  real  and  imaginary  parts. 

The  imaginary  part  of  the  equation  of  motion,  Equa¬ 
tion  (5),  essentially  determines  the  neutral  stability  condition 
of  the  system,  and  the  real  part  determines  the  frequency  of 
oscillation.  Of  course  now  both  of  these  results  depend  on  the 
pitch  amplitude  of  motion. 

While  structural  damping  is  readily  included  in  the  analy¬ 
sis,  as  will  be  seen  below,  it  will  be  helpful  to  understand  the 
essence  of  the  results  by  first  considering  the  solution  for  zero 
structural  damping. 

Zero  Structural  Damping 

In  this  case,  Equation  (5)  states  that  a  neutrally  stable  oscil¬ 
lation  will  occur  when  the  imaginary  part  of  the  aerodynamic 
moment  becomes  zero.  This  will  occur  at  some  reduced  fre¬ 
quency  for  a  particular  angle  of  pitch  oscillation  (and  other 
parameters  fixed  such  as  Mach  number).  Then  from  Equa¬ 
tion  (4),  one  can  solve  for  the  frequency  of  this  neutrally  sta¬ 
ble  oscillation.  For  sufficiently  small  motions,  this  is  the  flutter 
solution;  for  larger  motions,  we  determine  a  limit  cycle  oscil¬ 
lation.  The  solution  procedure  then  is  to  select  an  amplitude 
of  oscillation,  determine  the  reduced  frequency  at  which  the 
imaginary  part  of  the  aerodynamic  moment  is  zero  from  Fig¬ 
ure  (5),  and  then  determine  the  frequency  of  the  oscillation 
from  Equation  (4).  Note  this  is  essentially  the  same  computa¬ 
tional  procedure  as  for  a  classical  flutter  solution,  except  now 
the  reduced  frequency,  the  frequency  of  oscillation,  and  the 
reduced  velocity  are  all  functions  of  the  pitch  amplitude. 

It  should  be  noted  however  that  just  because  the  imaginary 
part  of  the  aerodynamic  moment  vanishes  (i.e.  the  aerody¬ 
namic  damping  becomes  zero),  that  alone  does  not  insure  that 
a  neutrally  stable  oscillation  will  occur.  This  is  because  the 
frequency  determined  from  Equation  (4)  must  be  physically 
possible,  i.e  the  right  hand  side  of  Equation  (4)  must  be  posi¬ 
tive.  It  is  evident  that  the  right  hand  side  of  Equation  (4)  de¬ 
pends  only  on  the  reduced  frequency  (which  is  known  by  the 
requirement  that  the  imaginary  part  of  the  aerodynamic  mo¬ 
ment  be  zero)  and  a  non-dimensional  moment  of  inertia.  But 
of  course  these  reduced  frequencies  themselves  depend  on  the 
pitch  amplitude.  Thus  one  can  determine  when  the  right  hand 
side  of  Equation  (4)  is  positive  or  negative  and  express  the  re¬ 
sult  in  terms  of  pitch  amplitude  and  moment  of  inertia.  This 
relationship  is  shown  in  Figure  (6),  and  the  regions  where  flut¬ 
ter  and  limit  cycle  oscillations  are  or  are  not  possible  are  indi¬ 
cated.  The  value  of  moment  of  inertia  that  marks  the  boundary 
between  no  flutter  or  LCO  possible  and  possible  flutter  or  LCO 
is  termed  the  ’'asymptotic  value”. 

Large  Pitch  Moment  of  Inertia 

Now  if  the  mass  ratio  or  moment  of  inertia  is  much  larger 
than  the  asymptotic  value,  a  not  uncommon  circumstance,  then 
the  flutter  or  LCO  frequency  is  simply  equal  to  the  structural 
pitch  natural  frequency.  See  Equation  (4).  With  this  approx¬ 
imation,  the  results  of  Figure  (7)  are  obtained  for  both  zero 
and  non-zero  structural  damping.  Note  that  the  curves  bend  to 
the  left  which  is  indicative  of  an  unstable  LCO.  That  is,  these 
results  are  to  be  interpreted  as  the  amplitude  of  a  disturbance 
required  to  initiate  explosive  flutter  below  the  classical  flutter 


Figure  6:  Asymptotic  Value  of  Pitch  Inertia  For  Various  Pitch 
Amplitudes  Marking  Regions  Where  Flutter  and  Limit  Cycle 
Oscillations  Are  or  Are  Not  Possible:  NACA  64A010A  Air¬ 
foil  Section,  M  =  0.80,  o0  =  0.0  (deg),  Elastic  Axis  Location, 
a  -  e/b  -  -0.6 
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Figure  7:  LCO  Amplitude  Versus  Reduced  Velocity:  NACA 
64A010A  Airfoil  Section,  M  =  0.80,  cx0  =  0.0  (deg),  Elastic  Axis 
Location,  a  =  e/b  =  -0.6 


velocity  for  this  single-degree-of-freedom  pitch  oscillation. 

In  Figure  (8),  shown  are  the  values  of  structural  damping 
(normalized  by  pitch  moment  of  inertia)  that  correspond  to 
neutrally  stable  limit  cycle  oscillations.  These  can  be  calcu¬ 
lated  from  Equation  (5)  as  a  function  of  reduced  velocity  for 
various  pitch  amplitudes.  A  cross-plot  of  these  data  is  used  to 
construct  the  plots  for  non-zero  damping  values  as  shown  in 
Figure  (7). 
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Figure  8:  Normalized  Structural  Damping  Corresponding  to 
Neutrally  Stable  Limit  Cycle  Oscillations  for  Large  Pitch  Mo¬ 
ment  of  Inertia:  NACA  64A010A  Airfoil  Section,  M  =  0.80, 
qo  =  0.0  (deg),  Elastic  Axis  Location,  a  =  e/b  =  -0.6 


Effects  of  Finite  Pitch  Moment  of  Inertia 

For  general  values  of  moment  of  inertia  and  structural 
damping,  the  solution  algorithm  using  Equations  (4)  and  (5) 
proceeds  as  follows.  First  select  a  Mach  number  and  pitch 
axis,  and  for  a  range  of  pitch  amplitudes,  determine  the  first 
harmonic  of  the  aerodynamic  moment  (including  higher  har¬ 
monics  of  the  aerodynamic  model  and  their  effect  on  the  fun¬ 
damental  harmonic).  Then  for  a  given  pitch  amplitude,  choose 
a  reduced  frequency  and  determine  the  flutter  or  LCO  oscil¬ 
lation  frequency  from  Equation  (4).  This  frequency  will  be 
proportional  to  the  pitch  structural  frequency,  of  course.  With 
the  flutter  or  LCO  frequency  determined,  and  the  reduced  fre¬ 
quency  selected,  one  then  knows  the  flow  velocity  correspond¬ 
ing  to  the  chosen  pitch  amplitude.  Finally  from  Equation  (5), 
determine  the  structural  damping  value  necessary  to  give  a 
neutrally  stable  flutter  or  limit  eyrie  oscillation.  From  this 
perspective,  the  flutter  condition  is  simply  the  neutrally  stable 
motion  that  may  exist  at  small  angles  of  twist,  and  the  LCO 
are  the  neutrally  stable  oscillations  that  may  exist  when  the 
pitch  amplitude  is  finite.  Of  course  the  flutter  or  LCO  may 
become  unstable  when  it  is  perturbed  (e.g.  by  perturbations  in 
the  amplitude  of  oscillation),  and  this  is  indeed  the  case  in  the 
example  treated  here. 

Up  to  this  point,  we  have  assumed  that  the  pitch  moment 
of  inertia  is  well  above  its  asymptotic  value.  Hence  the  flutter 
frequency  is  the  same  as  the  structural  natural  pitch  frequency. 

Now  we  consider  the  more  general  case  and  a  range  of 
pitch  inertias  such  that  the  flutter  frequency  is  no  longer 
precisely  equal  to  the  structural  natural  frequency  in  pitch. 
Results  are  shown  for  non-dimensional  pitch  inertias  of 
200,100,50,37.5  and  25  in  Figures  (9)  and  (10).  These  are  for 
LCO  amplitude  and  frequency,  respectively,  versus  reduced 
velocity.  The  asymptotic  pitch  inertia  results  are  also  shown 
for  reference. 

As  expected,  for  sufficiently  large  pitch  inertia,  say  greater 
than  200,  the  asymptotic  results  are  good  approximations. 


Figure  9:  LCO  Amplitude  Versus  Reduced  Velocity  for  Vari¬ 
ous  Pitch  Inertias:  NACA  64A010A  Airfoil  Section,  M  =  0.80, 
qo  =  0.0  (deg),  Elastic  Axis  Location,  a  =  e/b  =  -0.6 


Figure  10:  LCO  Amplitude  Versus  LCO  Frequency  for  Vari¬ 
ous  Pitch  Inertias:  NACA  64A010A  Airfoil  Section,  M  =  0.80, 
qo  a  o.O  (deg),  Elastic  Axis  Location,  a  =  e/b  =  -0.6 


However  for  pitch  inertias  less  than  100,  the  results  show  a 
more  sensitive  dependence  on  pitch  moment  of  inertia.  For 
sufficiently  small  pitch  moment  of  inertia,  of  course,  no  flutter 
or  LCO  is  possible. 

The  relationship  between  pitch  moment  of  inertia  and  re¬ 
duced  velocity  may  be  even  more  clearly  seen  by  fixing  the 
pitch  amplitude  and  then  plotting  these  variables  as  shown  in 
Figure  (1 1).  Note  in  this  figure,  as  reduced  velocity  decreases, 
the  pitch  moment  of  inertia  for  flutter  and  LCO  to  occur  tends 
to  infinity.  Thus  for  sufficiently  small  reduced  velocity  no  flut¬ 
ter  or  LCO  will  occur.  Conversely  as  the  pitch  moment  of  iner¬ 
tia  decreases,  the  reduced  velocity  for  flutter  or  LCO  to  occur 
tends  to  infinity.  Thus  below  some  value  of  pitch  moment  of 
inertia,  no  flutter  or  LCO  is  possible.  Of  course,  these  results 


6 


Figure  11:  Pitch  Inertia  Versus  Reduced  Velocity  for  Fixed 
Pitch  Amplitudes:  NACA  64A010A  Airfoil  Section,  M  =  0.80, 
Q0  =  0.0  (deg),  Elastic  Axis  Location,  a  =  e/b  =  -0.6 


are  for  a  fixed  pitch  amplitude  when  in  fact  the  pitch  ampli¬ 
tude  is  an  outcome  of  the  analysis  (not  an  input).  However, 
the  results  are  not  very  sensitive  to  pitch  amplitude  and  the 
conclusions  regarding  asymptotic  behavior  hold  over  the  full 
range  of  pitch  amplitudes  considered  here. 

CONCLUSIONS 

Aerodynamic  nonlinearities  may  be  give  rise  to  LCO, 
and  these  may  be  either  stable  (favorable)  or  unsta- 
ble(unfavorable).  An  example  of  the  former  is  shown  here 
as  the  nonlinear  counterpart  of  classical  linear  aeroelastic  di¬ 
vergence.  An  example  of  the  latter  is  also  shown  here  as  the 
nonlinear  counterpart  of  single-degree-of-ffeedom  pitch  flut¬ 
ter.  Future  work  will  be  directed  toward  the  study  of  the  non¬ 
linear  counterpart  of  classical  bending/torsion  flutter  where 
similar  methods  may  be  used  and  results  for  the  two-degree- 
of-freedom  case  will  be  presented  in  the  final  paper. 
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INTRODUCTION 

In  the  following  extended  abstract,  we  demonstrate  how 
the  recently  devised  proper  orthogonal  decomposition  (POD) 
based  reduced  order  modeling  (ROM)  technique  (Refer¬ 
ences  [1]  and  [2])  can  be  used  to  model  unsteady  aerodynamic 
and  aeroelastic  characteristics  of  three-dimensional  transonic 
wing  configurations.  Although  transonic  Euler  flows  are 
considered  in  [1]  and  [2],  the  initial  demonstrations  of  the 
POD/ROM  method  as  presented  in  these  references  are  for 
two-dimensional  flow  and  two  structural  degree-of-freedom 
airfoil  configurations.  Also  in  [3],  an  application  of  the 
POD/ROM  technique  to  the  well  known  vortex  lattice  method 
has  been  presented. 

In  extending  the  POD/ROM  technique  to  three- 
dimensions,  two  primary  issues  have  been  of  concern. 
First,  the  size  of  the  computational  fluid  dynamic  (CFD) 
model  will  in  general  be  at  Least  an  order  of  magnitude  greater 
than  for  two-dimensions.  Whereas  a  typical  CFD  model  for 
a  realistic  two-dimensional  configuration  might  have  on  the 
order  of  1 0’s  or  even  1 00's  of  thousands  of  degrees  of  freedom 
(DOF),  a  CFD  model  for  a  three-dimensional  configuration 
might  easily  have  on  the  order  of  at  least  hundreds  of  thou¬ 
sands  if  not  millions  or  more  DOF’s.  In  two-dimensions,  we 
have  found  that  very  accurate  ROM’s  with  on  the  order  of 
only  a  few  dozen  DOF’s  can  be  devised  using  the  POD/ROM 
methodology.  A  first  issue  to  address  has  thus  been  whether 
or  not  in  three-dimensions  one  can  also  generate  accurate 
ROM’s,  which  require  at  most  a  few  dozen  DOF’s. 

The  second  concern  is,  for  any  variation  of  the  structural 
properties  of  a  given  wing  under  consideration,  will  a  com¬ 
pletely  new  ensemble  of  solution  vector  “snapshots”  have  to 
be  computed  in  order  to  devise  an  accurate  POD/ROM.  A  ba¬ 
sic  aspect  of  the  POD/ROM  method  is  that  an  ensemble  of  so¬ 
lution  vectors  is  first  assembled  by  computing  unsteady  CFD 
solutions  at  a  number  of  discrete  frequencies  within  a  fre¬ 
quency  range  of  interest  for  the  unsteady  structural  motions 
that  are  also  of  interest.  In  two  dimensions,  this  step  is  rel¬ 
atively  straight  forward  since  one  only  has  to  consider  a  few 
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possible  motions,  e.g.  pitch  and  plunge. 

In  three-dimensions  however,  the  wing  vibratory  mode 
shapes  will  be  different  for  each  different  structural  configu¬ 
ration  of  a  given  wing.  There  can  in  fact  be  an  infinite  num¬ 
ber  of  unsteady  motions  (or  at  least  a  substantial  number  of 
motions  equivalent  to  the  number  of  DOF’s  ol  the  discrete 
structural  model).  Thus  the  second  concern  about  extending 
the  POD/ROM  to  three-dimensions  has  been  whether  or  not 
it  is  necessary  to  compute  a  completely  different  ensemble  of 
solution  snapshots  for  every  possible  structural  configuration. 
For  example,  say  one  computes  solution  snapshots  for  a  given 
wing  configuration  based  on  the  wing’s  particular  vibratory 
modes  shapes  in  order  to  develop  a  POD/ROM  to  model  the 
configuration’s  aeroelastic  characteristics.  Then  the  question 
is,  if  the  structural  make-up  of  the  wing  changes,  does  one 
have  to  compute  a  whole  new  ensemble  of  solution  snapshots 
for  the  same  wing,  but  for  the  different  set  of  vibratory  modes 
shapes. 

Fortunately  in  addressing  these  two  issues,  and  as  will 
be  shown  in  the  following,  we  have  found  that  accurate 
POD/ROM’s  with  just  a  few  dozen  degrees  of  freedom  can 
in  fact  be  created  for  a  realistic  transonic  three-dimensional 
configurations.  This  is  true  even  though  in  the  model  prob¬ 
lem  to  be  shown  subsequently,  the  CFD  model  is  easily  an 
order  of  magnitude  larger  than  anything  we  have  previously 
studied  in  two-dimensions.  Furthermore,  we  have  discovered 
that  a  “fundamental”  ensemble  of  solution  snapshots,  based  on 
wing  motions  that  need  not  be  related  to  the  structural  moues 
under  consideration,  can  be  assembled  as  a  first  step.  Accu¬ 
rate  POD/ROM’s  for  a  given  wing  configuration  can  then  be 
created  by  simply  adding  to  this  “fundamental”  ensemble,  the 
snapshots  corresponding  to  actual  wing  structural  modal  mo¬ 
tions  solely  at  the  frequencies  corresponding  to  the  end  points 
of  the  frequency  range  of  interest.  In  general,  these  two  snap¬ 
shots  prove  to  be  sufficient  to  “lock  in”  the  conditions  corre¬ 
sponding  to  the  particular  structural  motion,  and  indeed  the 
fundamental  ensemble  of  solution  snapshots  is  sufficient  to  re¬ 
veal  the  unsteady  dynamics  of  the  fluid  dynamic  model.  The 
fundamental  ensemble  of  snapshots  can  be  used  again  and 
again  even  as  the  structural  mode  change,  and  thus  the  com¬ 
putational  cost  of  having  to  compute  an  entirely  new  snapshot 
ensemble  for  every  new  structural  configuration  is  greatly  re¬ 
duced. 

POD/ROM  METHODOLOGY 

In  the  following,  we  will  be  considering  inviscid  three- 
dimensional  Euler  flows.  More  specifically,  linearized  (about 
some  nonlinear  background  steady  flow)  unsteady  frequency- 
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domain  CFD  solutions  to  the  Euler  equations  are  computed. 
The  POD/ROM  procedure  can  be  considered  as  a  “wrapper” 
around  any  typical  CFD  method,  and  the  CFD  method  we  have 
employed  for  the  present  analysis  is  a  variant  of  Ni’s  [4]  ap¬ 
proach  to  the  standard  Lax-Wendroff  method.  The  frequency 
domain  CFD  method  in  effect  represents  a  linear  system  for¬ 
mulation  of  the  unsteady  fluid  dynamic  model,  i.e. 

Aq=-BC  (1) 


Here,  #  is  an  N  x  K  matrix  whose  fcth  column  is  the  shape 
vector  <£fc,  and  £  is  the  K  dimensional  vector  of  augmented 
aerodynamic  state  variables 

A  reduced-order  representation  of  the  fluid  dynamic  and 
aeroelastic  systems  can  be  formulated  by  substituting  Equa¬ 
tion  (5)  into  Equation  (1)  and/or  (3)  and  pre-multiplying  by 
the  Hermitian  transpose  ($/f)  of  i.e. 

=  $hBC  or  At  =  -BC  (6) 


where  q  is  an  N  dimensional  vector  (N  is  the  number  of  mesh 
points  time  the  number  of  dependent  flow  variables)  of  the  un¬ 
known  flow  variables  a  each  mesh  point  in  the  CFD  domain, 
and  C  is  the  L  dimensional  vector  (L  is  the  number  structural 
mode  shapes)  of  modal  coordinates  for  the  structural  model. 
A  is  the  TV  x  N  fluid  dynamic  influence  matrix,  and  B  is  the 
N  x  L  matrix  which  relates  the  flow  solver  boundary  condi¬ 
tions  to  each  particular  mode  shape.  Both  A  and  B  are  func¬ 
tions  of  the  background  flow  and  unsteady  frequency  cj.  The 
structural  equations  for  the  wing  configuration  being  modeled 
within  the  flow  can  be  written  as 

DC  =  -Cq  (2) 

where  D  represents  the  L  x  L  structural  influence  matrix  (i.e. 
D  =  - u)2M,  +  K  where  M  and  K  are  the  generalized  mass 
and  stiffness  matrices),  and  C  is  the  L  x  M  matrix  which 
represents  the  discrete  integration  used  to  obtain  the  general¬ 
ized  forces  associated  with  each  modes  shape  based  on  the 
unsteady  flow  q.  When  Equations  (1)  and  (2)  are  put  together, 


The  resulting  equation  (3)  is  a  fully  coupled  aeroelastic  sys¬ 
tem  of  equations,  which  for  nontrivial  q  and  C,  represents  an 
eigenvalue  problem  with  u>  being  the  eigenvalue.  Any  eigen¬ 
values  with  a  positive  real  part  imply  the  aeroelastic  system  is 
unstable. 

The  problem  with  constructing  and  solving  this  eigenvalue 
problem  is  that  A  is  simply  too  large  for  realistic  configura¬ 
tions.  As  mentioned  in  the  introduction,  N  can  easily  be  on 
the  ji  der  of  10,000  or  100,000  for  two-dimensional  configura¬ 
tions,  and  on  order  of  100,000  to  1,000,000  or  even  more  for 
three  dimensional  configurations.  For  such  large  cases,  even 
attempting  to  set  up  A  is  well  beyond  the  memory  limits  of 
todays  largest  computers. 

The  basic  premise  of  the  POD/ROM  methodology  is  that 
we  assume  the  unknown  flowfield  solution  vector  q  can  be 
expressed  as  a  Ritz  type  expansion  of  the  form 

*«JV  (4) 


where  tk  is  a  generalized  coordinate  sometimes  referred  to  as 
an  augmented  aerodynamic  state  variable,  and  (p^  is  the  cor¬ 
responding  Ritz  vector.  Equation  (4)  can  also  be  written  in 
matrix  form  as 


and 


If  the  Ritz  approximation  is  good  one,  ( K  <$C  N)y  and  Equa¬ 
tions  (6)  and  (7)  will  represent  much  smaller  systems  that  can 
readily  be  solved  using  conventional  eigenvalue  techniques. 

The  next  question  becomes  what  are  good  choices  for  the 
Ritz  vectors  4>k  that  will  in  fact  result  in  good  Ritz  approxi¬ 
mations.  Previous  studies  as  detailed  in  References  [1]  and  [2] 
have  demonstrated  that  shape  vectors  derived  via  the  proper 
orthogonal  decomposition  technique  (see  for  instance  [5],  [6], 
and  [7])  are  an  excellent  source.  For  the  sake  of  brevity,  the 
of  details  are  omitted  here,  but  a  discussion  of  how  the  shapes 
are  derived  can  be  found  in  References  [1]  and  [2].  The  ba¬ 
sic  premise  behind  their  formulation  is  that  a  number  solu¬ 
tion  “snapshots”  are  directly  computed  for  a  number  of  dis¬ 
crete  frequencies  and  unsteady  structural  motions  of  interest. 
From  this  ensemble  of  solution  vectors,  the  POD  shapes  are 
easily  derived  by  solving  a  small  (the  size  of  the  number  of 
snapshots)  eigenvalue  problem.  The  first  few  POD  modes  de¬ 
scribe  the  most  dominant  dynamic  characteristics  of  the  fluid 
dynamic  system,  and  as  such,  the  POD  shapes  have  proven 
to  be  an  excellent  set  of  Ritz  vectors  for  fluid  dynamic  and/or 
aeroelastic  models. 

MODEL  PROBLEM 

The  configuration  under  consideration  is  the  AGARD 
model  445.6  wing  (Reference  [8]  and  [9]).  This  is  a  45  degree 
quarter  chord  swept  wing  using  the  NACA  64A004  airfoil  sec¬ 
tion  that  has  an  aspect  ratio  of  3.3  (for  the  full  span)  and  a  taper 
ratio  of  2/3.  Figure  (1 )  illustrates  the  computational  mesh  em¬ 
ployed  for  this  configuration.  The  grid  is  an  “O-O”  topology 
that  employs  49  computational  nodes  about  each  airfoil  sec¬ 
tion,  33  nodes  normal  to  the  wing,  and  33  nodes  along  the 
semispan.  The  outer  boundary  of  the  grid  extends  five  semis¬ 
pans  from  the  midchord  of  the  wing  root  section.  The  particu¬ 
lar  structural  configuration  of  the  wing  is  referred  to  as  the  2.5 
ft.  weakened  model  3  (again  see  Reference  [8]  and  [9]). 

Figure  (2)  shows  the  computed  wing  surface  and  symme¬ 
try  plane  steady  flow  pressure  contours  for  the  Mach  numbers 
of  0.960  and  1.141.  A  relatively  weak  shock  can  be  seen 
at  the  trailing  edge  for  M=0.960.  This  shock  appears  to  get 
stronger  at  M=1.141.  The  wing  section  is  quite  thin  (4%),  so  a 
strong  shock  is  not  really  expected.  Comparing  contours,  our 
flowfields  look  very  comparable  to  those  of  Lee-Rausch  and 
Batina  [10],  although  they  employed  a  much  larger  mesh. 

FLUTTER  RESULTS 

Figure  (3)  shows  the  eigenvalue  root-loci  when  sweeping 
through  various  mass  ratios  (to  which  there  is  a  correspond- 
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(b)  Outer  Boundary  Grid 


(b)  M=1.141 


Figure  1 :  AGRAD  445.6  W'lng  Grid  Topology 


ing  flutter  speed  index)  when  solving  the  aeroelastic  eigen¬ 
value  problem  posed  by  the  reduced-order  aeroelastic  model 
(Equation  (7)  for  various  Mach  numbers.  Solution  snapshots 
are  computed  for  the  first  five  given  wing  mode  shapes  for  re¬ 
duced  frequencies  (k  =  tob/Uoo)  from  k  =  0.0  to  k  =  0.5  in 
A k  =  0.1  increments.  This  configuration  flutters  for  frequen¬ 
cies  less  than  0.5,  and  as  such,  solution  snapshots  for  k  >  0.5 
are  unnecessary.  This  results  in  a  total  of  55  available  POD 
shape  vectors.  In  Figure  (3),  the  curves  represent  the  eigen¬ 
values  corresponding  to  the  primarily  structural  natural  modes 
as  mass  ratio  is  varied.  Our  method  also  determines  the  aeroe¬ 
lastic  modes  originating  from  the  fluid  dynamic  modes  of  the 
POD/ROM.  For  the  range  of  mass  ratios  (0  <  /z  <  500)  swept 
through  in  these  parametric  analyses,  the  fluid  dynamic  modes 
are  very  damped,  and  as  such  lie  to  the  left  and  outside  of  the 
eigenspectrum  range  we  show.  As  can  been,  for  each  of  the 
Mach  numbers,  the  first  structural  mode  tends  to  be  the  criti- 


Flgure  2:  AGRAD  445.6  Wing  Surface  and  Symmetry  Plane 
Pressure  Contours 

cal  flutter  mode.  For  the  highest  Mach  number  however,  the 
third  structural  mode  can  go  unstable  if  the  mass  ratio  is  large 
enough.  Also  from  this  figure,  it  is  evident  that  it  is  unneces¬ 
sary  to  use  all  55  of  the  available  POD  shapes.  If  fact,  with  less 
than  one  half  of  the  POD  modes  (25  for  instance),  relatively 
converged  results  (in  the  sense  of  POD  mode  refinement)  can 
be  achieved. 

Figure  (4)  shows  the  computed  POD/ROM  flutter  speed 
and  flutter  frequency  ratios,  along  with  experimental  data  [8], 
and  data  from  two  other  computational  methods  ([10] 
and  [11])  ,  as  a  function  of  Mach  number.  As  can  be  seen, 
using  our  methodology,  we  produce  the  well  known  transonic 
flutter  speed  dip,  and  our  results  are  all  within  the  same  tol¬ 
erance  to  the  experimental  results  as  the  other  computational 
methods.  Gupta  [11]  does  show  better  agreement  with  ex¬ 
periment  at  the  two  supersonic  Mach  numbers,  and  Gupta  at¬ 
tributes  this  better  agreement  to  better  CFD  grid  refinement. 
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Figure  3:  Aeroelastic  Root-Loci  at  Various  Mach  Numbers 
for  the  AGRAD  445.6  Wing  “Weakened”  Configuration,  a0 
=  0.0  (deg) 


In  the  final  paper,  we  will  also  address  this  issue. 

Figure  (5)  again  shows  the  computed  POM/ROM  flutter 
speed  and  flutter  frequency  ratios  as  a  function  of  Mach  num¬ 
ber.  In  this  instance  however,  the  flutter  results  are  shown  for 
various  fractions  of  the  total  available  POD  shapes.  As  can 
be  seen,  even  with  as  few  as  1/2  of  the  available  shapes,  rela¬ 
tively  well  converged  results  are  obtained.  Note  however  there 
is  somewhat  greater  sensitivity  to  the  number  of  POD/ROM 
modes  retained  at  the  supersonic  Mach  numbers. 

THE  USE  OF  ALTERNATE  MODAL 
EXCITATIONS  FOR  SNAPSHOTS 

One  of  the  key  concerns  towards  in  the  POD/ROM  method 
to  three-dimensional  flows  has  been  whether  or  not  an  entire 
set  of  solution  snapshots  must  be  computed  for  each  possible 
structural  configuration  of  interest.  That  is,  say  we  wish  to 
consider  a  similarly  shaped  wing  that  has  a  slightly  different 


Mach  Number,  M 

(a)  Flutter  Speed 


Mach  Number,  M 


(b)  Flutter  Frequency  Ratio 


Figure  4:  Mach  Number  Flutter  Trend  for  the  AGRAD  445.6 
Wing  “Weakened"  Configuration,  a0  =  0.0  (deg) 


structural  configuration,  which  in  turn  means  the  wing  vibra¬ 
tory  mode  shapes  are  different.  Does  this  mean  that  one  has 
to  go  through  and  compute  a  whole  new  ensemble  of  solu¬ 
tion  snapshots  based  on  these  new  structural  motions  in  order 
to  do  a  flutter  analysis  for  the  new  wing  configuration.  For¬ 
tunately,  as  we  will  demonstrate  in  the  following,  although  a 
few  snapshots  based  on  the  new  modal  motion  will  need  to 
be  computed,  the  larger  number  of  snapshots  computed  at  nu¬ 
merous  frequencies  will  be  unnecessary.  The  snapshots  com¬ 
puted  from  a  previous  wing  structural  configuration  will  still 
serve  the  purpose.  That  is,  a  couple  of  solution  snapshots  will 
be  needed  for  the  new  structural  motions,  however,  these  will 
only  need  to  be  computed  at  the  end  points  of  the  frequency 
range  of  interest.  Figure  (6)  demonstrates  how  this  works. 

Figure  (6a)  shows  the  real  and  imaginary  parts  of  the  coef- 
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(b)  Flutter  Frequency  Ratio 

Figure  5:  POD/ROM  Shape  Vector  Refinements  Characteristic 
for  Mach  Number  Flutter  Trend  for  AGRAD  445.6  Wing  “Weak¬ 
ened"  Configuration,  ao  =  0.0  (deg) 


flcient  of  the  generalized  aerodynamic  force  corresponding  to 
the  first  mode  pressure  acting  through  the  first  mode  shape  as  a 
function  of  reduced  frequency  at  a  Mach  number  of  0.960.  The 
coefficient  of  the  generalized  aerodynamic  force  is  defined  as 

(k)  — - 2  [  f  Pj(k)nz  dA  (8) 

Qoo^r  J  J 

where  =  Poo^/2  is  the  freestream  dynamic  pressure 
and  Cr  is  the  root  chord  length.  The  integral  is  evaluated  over 
the  surface  of  the  wing,  and  nz  is  z  component  of  the  wing 
surface  normal  vector  (i.e.  n  =  nx i  T*  ftyj  +  ^zk  and  n  is 
oriented  to  point  towards  the  wing  surface).  In  this  definition, 
Pj(k)  represents  the  frequency  dependent  unsteady  pressure 


resulting  from  a  wing  deformation  motion  of 

-  -  4>i  (9) 

Cr 

The  curves  presented  in  Figure  (6a)  are  based  on  the  actual 
solution  snapshots  and  thus  are  what  we  desire  the  POD/ROM 
to  model.  In  Figure  (6b),  the  POD/ROM  of  CQl  l  based  on 
snapshots  for  each  ot  the  five  structural  mode  shapes  at  fre¬ 
quencies  of  A;  =  0.0  and  fc  =  1.0  (for  a  total  of  10  snapshots) 
is  compared  against  Cqx  ,  for  the  actual  snapshots  for  all  fre¬ 
quencies  between  k  =  0.0  and  k  —  1.0.  As  can  be  seen,  the 
POD/ROM  matches  at  the  end  points  of  the  frequency  range 
as  is  expected,  however  this  ca*de  POD/ROM  performs  rather 
poorly  for  the  intermediate  frequencies.  Of  course,  if  we  use 
snapshots  at  all  the  frequencies  between  k  =  0.0  and  k  =  1.0, 
the  POD/ROM  would  exactly  reproduce  the  data. 

Next,  in  Figure  (6c),  a  new  POD/ROM  for  CQx  i  now 
based  on  solution  snapshots  unrelated  to  the  actual  mode 
shapes  is  shown.  The  simple  wing  motion  snapshots  are  for  a 
full  wing  plunge  motion  (up/down),  full  wing  pitch  about  the 
quarter  chord,  a  first  bending  type  of  motion  (wing  is  fixed  at 
the  root,  and  the  z  coordinate  component  of  deflection  varies 
linearly  with  span),  and  a  first  twist  type  of  motion  (wing  is 
fixed  at  the  root,  and  the  pitch  varies  linearly  with  span)  for 
frequencies  from  k  =  0.0  to  k  =  1.0  at  A  A;  =  0.1  incre¬ 
ments  for  a  total  of  44  solution  snapshots.  As  can  be  seen, 
the  POD/ROM  in  this  case  also  perform  very  poorly.  Unbe¬ 
knownst  however,  these  solutions  are  in  fact  helping  to  reveal 
the  dynamics  of  the  system.  In  fact,  when  one  uses  these  snap¬ 
shots  in  combination  with  the  actual  structural  mode  snapshots 
solely  at  the  end  points  of  the  frequency  range  of  interest, 
one  gets  a  POD/ROM  which  produces  very  accurate  results 
to  Cqx1  as  is  evident  from  Figure  (6d). 

Figure  (7)  shows  a  comparison  of  the  POD/ROM  Mach 
number  flutter  trends  in  the  case  where  first,  the  POD/ROM  is 
based  on  solution  snapshots  corresponding  to  the  actual  modal 
shapes  of  the  wing  to  the  case  where  second,  the  POD/ROM  is 
based  on  snapshots  using  the  simple  wing  motion  mode  shapes 
as  discussed  in  the  previous  paragraph.  As  for  the  previous 
Mach  number  flutter  results  (Figures  (4)  and  (5),  the  snapshot 
reduced  frequencies  range  from  k  =  0.0  to  k  =  0.5  in  A  A;  = 
0.1  increments.  Included  In  this  second  ensemble  of  solution 
snapshots,  are  the  snapshots  corresponding  to  the  particular 
motions  of  the  actual  mode  shapes  at  the  end  points  of  the 
frequency  range  of  interest. 

As  can  be  seen  in  Figure  (7),  one  can  obtain  accurate 
POD/ROM  flutter  results  using  solution  snapshots  unrelated 
to  the  actual  wing  motions  (except  at  the  end  points  of  the  fre¬ 
quency  range  of  interest)  that  compare  very  well  to  the  flutter 
results  based  on  a  POD/ROM  model  using  snapshots  corre¬ 
sponding  to  the  actual  motions.  This  is  especially  true  at  the 
lower  Mach  numbers.  There  is  some  difference  at  the  high¬ 
est  Mach  number,  again  suggesting  supersonic  flow  is  more 
sensitive  for  this  wing, 

To  further  illustrate  the  concept  of  being  able  to  use  solely 
the  frequency  range  of  interest  end  point  snapshots,  we  present 
a  few  results  for  a  simple  two-dimensional  configuration  that 
first  led  us  to  the  idea  of  the  using  the  technique  in  three- 
dimensions.  In  an  analogous  situation  of  having  to  consider 
multiple  structural  degrees  of  freedom  in  three-dimensions,  we 
first  considered  an  unsteady  NACA  64A010A  airfoil  configu- 
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ration  that  not  only  under  goes  typical  plunge  and  pitch  mo¬ 
tions  (see  Figures  (8a)  and  (8b),  but  also  has  motions  where 
airfoil  mean camber  line  distorts_based  on  simple  trigono¬ 
metric  functions,  i.e.  zc(x)  =  6X  cos(2t rx/c]  (Figure  8c), 
zc{x)  =  h  sin(27r^/c)  (Figure  8d),  zc(x)  =  S3  cos(4ttx/c) 
(Figure  8e),  etc.  The  initial  question  was,  as  one  considers 
each  subsequent  motion,  does  one  have  to  include  a  number  of 
snapshots  based  on  the  new  motion  that  is  equal  to  the  num¬ 
ber  of  snapshots  for  each  of  the  previous  motions  in  order  to 
produce  an  accurate  POD/ROM. 

Figure  (9)  illustrates  how  after  a  sufficient  number  of  snap¬ 
shots  have  been  included  in  the  snapshot  ensemble,  only  the 
end  point  frequencies  are  required  for  each  additional  motion. 
In  this  instance,  the  NACA  64A010A  airfoil  is  modeled  in  a 
M=0.5,qo  =  0.0  (deg)  background  flow,  and  shown  on  the  ab¬ 
scissa  are  increments  in  the  number  of  overall  motions  consid¬ 
ered.  Shown  on  the  ordinate  is  the  order  in  which  the  snapshots 
for  each  particular  motion  are  added  to  the  overall  ensemble. 
The  reduced  frequency  range  of  interest  is  0.0  <  k  <  1.0 
(note,  k  for  this  airfoil  problem  is  defined  as  k  =  cjc/E/qo), 
and  thus  the  first  two  snapshots  considered  for  each  motion 
correspond  the  end  points  of  this  frequency  range.  Further 
snapshots  added  to  the  ensemble  for  a  given  motion  are  done 
so  in  a  divide  and  conquer  strategy  to  best  model  the  interme¬ 
diate  frequencies. 

Considered  in  Figure  (9)  is  the  accuracy  of  modeling  the 
airfoil  unsteady  lift  and  moment  along  the  paths  s  =  rej9 
(where  0  =  90°,60o,  120°  and  0  <  r  <  1)  in  the  complex 
reduced  frequency  s  plane.  The  curves  illustrate  the  number 
of  snapshots  necessary  to  achieve  a  given  level  of  accuracy  for 
the  nth  and  all  previous  motions.  The  accuracy  is  based  on  a 
comparison  to  a  POD/ROM  that  is  derived  from  a  snapshot  en¬ 
semble  comprised  of  all  the  possible  motions  at  all  the  possible 
frequencies.  So  for  example,  to  achieve  a  10~3  L2  norm  ac¬ 
curacy  when  just  considering  plunge  motion,  one  needs  a  total 
of  ten  plunge  snapshots  of  the  frequency  values  indicated  on 
the  ordinate  of  the  plot.  If  next  considering  pitch  motion,  one 
would  then  need  to  add  only  three  pitch  motion  snapshots  cor¬ 
responding  to  the  frequencies  k  =  0.0,  k  =  1.0,  and  k  =  0.5 
to  the  overall  ensemble  to  get  10~3  L2  norm  accuracy  for  now 
both  the  pitch  and  plunge  motions.  Considering  next  the  first 
airfoil  bending  motion  S\ ,  one  would  then  need  to  add  a  total 
of  seven  snapshots  to  achieve  10“3  T/2_norm  accuracy  for 
now  plunge,  pitch,  and  motion.  Three  82  snapshots  would 
then  be  needed  when  also  taking  in  account  62  motions,  two 
83  snapshots  when  considering  83  motions,  and  so  on. 

As  can  be  seen,  Figure  (9)  illustrates  the  interesting  re¬ 
sult  that  after  a  sufficient  number  of  snapshots  have  been  in¬ 
cluded  in  the  overall  ensemble,  only  the  two  end  point  fre¬ 
quency  snapshots  for  each  subsequent  possible  motion  need 
be  added  to  the  ensemble  to  maintain  a  given  level  of  accu¬ 
racy.  Interestingly  enough,  this  appears  to  be  an  asymptotic 
limit.  That  is,  the  two  end  point  frequency  range  snapshots  al¬ 
ways  appear  to  be  necessary  when  considering  a  large  number 
of  possible  motions. 

CONCLUSIONS 

The  POD/ROM  method  has  been  demonstrated  for  the  flut¬ 
ter  analysis  of  a  three-dimensional  transonic  wing  configura¬ 
tion.  We  have  shown  that  the  number  of  ROM  DOF’s  neces¬ 


sary  to  create  accurate  models  is  on  the  order  of  a  few  dozen 
as  is  the  case  in  two-dimensions.  We  have  also  shown  that  it 
is  unnecessary'  to  compute  a  completely  new  ensemble  of  so¬ 
lution  snapshots  based  on  the  vibratory  mode  shapes  for  each 
new  structural  configuration  that  might  be  under  considera¬ 
tion.  One  can  simply  compute  a  set  of  snapshots  based  on 
some  basic  wing  motions  at  a  number  of  frequencies.  Then 
snapshots  only  at  the  end  points  of  the  frequency  range  of  in¬ 
terest  need  to  be  computed  for  the  specific  mode  shapes  of 
the  configuration  of  interest.  These  end  point  snapshots  “lock 
in”  the  unsteady  fluid  dynamic  characteristics  for  the  particu¬ 
lar  mode  shapes,  and  the  simple  motion  snapshots  then  act  to 
resolve  the  dominant  dynamics  of  the  flow  throughout  the  full 
frequency  range  of  interest. 

For  the  final  version  of  the  conference  paper,  will  we  in¬ 
clude  more  details  of  the  methodology  along  with  further  stud¬ 
ies  of  POD/ROM  refinements  including  mesh  convergence  is¬ 
sues. 
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(a)  Snapshots 


Reduced  Frequency,  k  =  ©b/U„ 

(b)  POD/ROM  Based  on  k  =  0.0  and  k  =  1.0 
Modal  Snapshots 


Reduced  Frequency,  k  =  cob/U_ 

(c)  POD/ROM  Based  on  “Basic"  Motions 
Snapshots 


Reduced  Frequency,  k  =  cob/U.. 

(d)  POD/ROM  Based  on  k  =  0.0  and  k  =  1.0 
Modal  Snapshots  and  “Basic”  Motions  Snap¬ 
shots 


(b)  Flutter  Frequency  Ratio 


Figure  7:  Mach  Number  Flutter  Trends  Using  Alternate  Unre¬ 
lated  Snapshots:  AGRAD  445.6  Wing  “Weakened”  Configura¬ 
tion,  oro  =  0.0  (deg) 


Figure  6:  Generalized  Force  Modeling  with  Unrelated  Mode 
Shape  Snapshots 
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